diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 5582fc70..bed7223d 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -74,7 +74,7 @@ jobs: strategy: matrix: os: [ "windows-latest", "ubuntu-latest" , "macos-latest"] - python-version: [ "3.9", "3.10", "3.11", "3.12" ] + python-version: [ "3.10", "3.11", "3.12", "3.13" ] runs-on: ${{matrix.os}} steps: - uses: actions/checkout@v4 diff --git a/.github/workflows/test_example.yml b/.github/workflows/test_example.yml index 0d4d1d96..877a49c6 100644 --- a/.github/workflows/test_example.yml +++ b/.github/workflows/test_example.yml @@ -7,7 +7,7 @@ jobs: strategy: matrix: os: [ "windows-latest", "ubuntu-latest" , "macos-latest"] - python-version: [ "3.9", "3.10", "3.11", "3.12"] + python-version: [ "3.10", "3.11", "3.12", "3.13" ] runs-on: ${{matrix.os}} steps: - uses: actions/checkout@v4 diff --git a/.github/workflows/test_pages.yml b/.github/workflows/test_pages.yml index ef949fc7..d4d33b9d 100644 --- a/.github/workflows/test_pages.yml +++ b/.github/workflows/test_pages.yml @@ -10,7 +10,7 @@ jobs: fetch-depth: 0 # otherwise, you will fail to push refs to dest repo - uses: actions/setup-python@v5 with: - python-version: '3.9' + python-version: '3.10' cache: 'pip' - name: Build and Commit uses: waltsims/pages@pyproject.toml-support diff --git a/.readthedocs.yaml b/.readthedocs.yaml index f6578e09..7837e0c0 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -3,7 +3,7 @@ version: 2 build: os: ubuntu-22.04 tools: - python: "3.9" + python: "3.11" sphinx: configuration: docs/conf.py python: diff --git a/examples/README.md b/examples/README.md index 16fada4d..bf5da62e 100644 --- a/examples/README.md +++ b/examples/README.md @@ -9,15 +9,15 @@ Every example has a short readme.md file which briefly describes the purpose of - [Array as a sensor](at_array_as_sensor/) ([original example](http://www.k-wave.org/documentation/example_at_array_as_sensor.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/at_array_as_sensor/at_array_as_sensor.ipynb)) - [Array as a source](at_array_as_source/) ([original example](http://www.k-wave.org/documentation/example_at_array_as_source.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/at_array_as_source/at_array_as_source.ipynb)) - [Linear array transducer](at_linear_array_transducer/) -([original example](http://www.k-wave.org/documentation/example_at_linear_array_transducer.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/at_linear_array_transducer/at_linear_array_transducer.ipynb)) + ([original example](http://www.k-wave.org/documentation/example_at_linear_array_transducer.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/at_linear_array_transducer/at_linear_array_transducer.ipynb)) - [Photoacoustic Waveforms](ivp_photoacoustic_waveforms/) ([original example](http://www.k-wave.org/documentation/example_ivp_photoacoustic_waveforms.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/ivp_photoacoustic_waveforms/ivp_photoacoustic_waveforms.ipynb)) - [Controlling the PML](na_controlling_the_pml/) -([original example](http://www.k-wave.org/documentation/example_na_controlling_the_pml.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/na_controlling_the_pml/na_controlling_the_pml.ipynb)) + ([original example](http://www.k-wave.org/documentation/example_na_controlling_the_pml.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/na_controlling_the_pml/na_controlling_the_pml.ipynb)) - [Defining An Ultrasound Transducer Example](us_defining_transducer) ([original example](http://www.k-wave.org/documentation/example_us_defining_transducer.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/us_defining_transducer/us_defining_transducer.ipynb)) - [Simulating Ultrasound Beam Patterns](us_beam_patterns/)([original example](http://www.k-wave.org/documentation/example_us_beam_patterns), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/us_beam_patterns/us_beam_patterns.ipynb)) - [Linear transducer B-mode](us_bmode_linear_transducer/) ([original example](http://www.k-wave.org/documentation/example_us_bmode_linear_transducer.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/us_bmode_linear_transducer/us_bmode_linear_transducer.ipynb)) - [Phased array B-mode](us_bmode_phased_array/) -([original example](http://www.k-wave.org/documentation/example_us_bmode_phased_array.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/us_bmode_phased_array/us_bmode_phased_array.ipynb)) + ([original example](http://www.k-wave.org/documentation/example_us_bmode_phased_array.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/us_bmode_phased_array/us_bmode_phased_array.ipynb)) - [Circular piston transducer](at_circular_piston_3D/) ([original example](http://www.k-wave.org/documentation/example_at_piston_and_bowl_transducers.php#heading3), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/at_circular_piston_3D/at_circular_piston_3D.ipynb)) - [Circular piston transducer (axisymmetric)](at_circular_piston_AS/) ([original example](http://www.k-wave.org/documentation/example_at_piston_and_bowl_transducers.php#heading4), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/at_circular_piston_AS/at_circular_piston_AS.ipynb)) @@ -32,6 +32,7 @@ Every example has a short readme.md file which briefly describes the purpose of - [Focussed Detector In 3D Example](sd_focussed_detector_3D/) ([original example](http://www.k-wave.org/documentation/example_sd_focussed_detector_3D.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/sd_focussed_detector_3D/sd_focussed_detector_3D.ipynb)) +- [Modelling Sensor Directivity In 2D Example](sd_directivity_modelling_2D/) ([original example](http://www.k-wave.org/documentation/example_sd_directivity_modelling_2D.php), [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/sd_directivity_modelling_2D/sd_directivity_modelling_2D.ipynb)) ## Contributing new examples @@ -39,13 +40,13 @@ When adding a new example notebook, follow these steps: 1. Search the [list of examples](https://docs.google.com/spreadsheets/d/1-x13iIez84AEyjjHMOe2GoC8FSyzUFHoy9R7VKttTlo/edit?usp=sharing) to find an open example. 1. Claim your example by opening an "Example" issue on GitHub. -3. Fork and clone the repository and create a branch for your new example. -2. Create an example sub-directory using the name from the hyperlink of the original k-wave example if it exists (e.g. for http://www.k-wave.org/documentation/example_ivp_loading_external_image.php name the directory "ivp_loading_external_image). -3. Add your example notebook to your example directory. -4. Create a Python script that mirrors your example notebook in the same directory. Using `nbconvert` is an easy way to convert to a Python script using ` jupyter nbconvert --to python --RegexRemovePreprocessor.patterns="^%" ` -5. Add a README.md file to your example directory briefly describing the concept or principle the example is meant to display and linking to the origonal k-wave example page if it exists. -6. Include a link in the `README.md` in the examples directory to a colab notebook for your example. -7. Add a your example to the list on this `README.md` and add a colab badge [using html](https://openincolab.com/) OR copy the pure markdown version above. -8. Open a pull request that [closes the open issue](https://docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue) from your forked example branch and name pull request "[Example] \". +1. Fork and clone the repository and create a branch for your new example. +1. Create an example sub-directory using the name from the hyperlink of the original k-wave example if it exists (e.g. for http://www.k-wave.org/documentation/example_ivp_loading_external_image.php name the directory "ivp_loading_external_image). +1. Add your example notebook to your example directory. +1. Create a Python script that mirrors your example notebook in the same directory. Using `nbconvert` is an easy way to convert to a Python script using ` jupyter nbconvert --to python --RegexRemovePreprocessor.patterns="^%" ` +1. Add a README.md file to your example directory briefly describing the concept or principle the example is meant to display and linking to the origonal k-wave example page if it exists. +1. Include a link in the `README.md` in the examples directory to a colab notebook for your example. +1. Add a your example to the list on this `README.md` and add a colab badge [using html](https://openincolab.com/) OR copy the pure markdown version above. +1. Open a pull request that [closes the open issue](https://docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue) from your forked example branch and name pull request "[Example] \". Thanks for contributing to k-wave-python! diff --git a/examples/sd_directivity_modelling_2D/README.md b/examples/sd_directivity_modelling_2D/README.md new file mode 100644 index 00000000..76e8c9d7 --- /dev/null +++ b/examples/sd_directivity_modelling_2D/README.md @@ -0,0 +1,7 @@ +# Modelling Sensor Directivity In 2D Example + +[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/waltsims/k-wave-python/blob/master/examples/sd_directivity_modelling_2D/sd_directivity_modelling_2D.ipynb) + +This example demonstrates how the sensitivity of a large single element detector varies with the angular position of a point-like source. + +To read more, visit the [original example page](http://www.k-wave.org/documentation/example_sd_directivity_modelling_2D.php). diff --git a/examples/sd_directivity_modelling_2D/sd_directivity_modelling_2D.ipynb b/examples/sd_directivity_modelling_2D/sd_directivity_modelling_2D.ipynb new file mode 100644 index 00000000..0c2848bc --- /dev/null +++ b/examples/sd_directivity_modelling_2D/sd_directivity_modelling_2D.ipynb @@ -0,0 +1,244 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install git+https://github.com/waltsims/k-wave-python " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Modelling Sensor Directivity In 2D Example\n", + "\n", + "This example demonstrates how the sensitivity of a large single element detector varies with the angular position of a point-like source." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "import os\n", + "from copy import deepcopy\n", + "from tempfile import gettempdir\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from kwave.data import Vector\n", + "from kwave.kgrid import kWaveGrid\n", + "from kwave.kmedium import kWaveMedium\n", + "from kwave.ksensor import kSensor\n", + "from kwave.ksource import kSource\n", + "from kwave.kspaceFirstOrder2D import kspaceFirstOrder2DC\n", + "from kwave.options.simulation_execution_options import SimulationExecutionOptions\n", + "from kwave.options.simulation_options import SimulationOptions\n", + "from kwave.utils.conversion import cart2grid\n", + "from kwave.utils.filters import filter_time_series\n", + "from kwave.utils.mapgen import make_cart_circle\n", + "from kwave.utils.matlab import ind2sub, matlab_find, unflatten_matlab_mask\n", + "from kwave.utils.colormap import get_color_map\n", + "from kwave.utils.data import scale_SI" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Create the computational grid and medium" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "grid_size_points = Vector([128, 128]) # [grid points]\n", + "grid_size_meters = Vector([50e-3, 50e-3]) # [m]\n", + "grid_spacing_meters = grid_size_meters / grid_size_points # [m]\n", + "kgrid = kWaveGrid(grid_size_points, grid_spacing_meters)\n", + "\n", + "# define the properties of the propagation medium\n", + "medium = kWaveMedium(sound_speed=1500)\n", + "\n", + "# define the array of time points [s]\n", + "Nt = 350\n", + "dt = 7e-8 # [s]\n", + "kgrid.setTime(Nt, dt)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define a large area detector" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sz = 20 # [grid points]\n", + "sensor_mask = np.zeros(grid_size_points)\n", + "sensor_mask[grid_size_points.x // 2, (grid_size_points.y // 2 - sz // 2) : (grid_size_points.y // 2 + sz // 2) + 1] = 1\n", + "sensor = kSensor(sensor_mask)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define a time varying sinusoidal source" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# define equally spaced point sources lying on a circle centred at the\n", + "# centre of the detector face\n", + "radius = 30 # [grid points]\n", + "points = 11\n", + "circle = make_cart_circle(radius * grid_spacing_meters.x, points, Vector([0, 0]), np.pi)\n", + "\n", + "# find the binary sensor mask most closely corresponding to the Cartesian\n", + "# coordinates from makeCartCircle\n", + "circle, _, _ = cart2grid(kgrid, circle)\n", + "\n", + "# find the indices of the sources in the binary source mask\n", + "source_positions = matlab_find(circle, val=1, mode=\"eq\")\n", + "\n", + "source = kSource()\n", + "source_freq = 0.25e6 # [Hz]\n", + "source_mag = 1 # [Pa]\n", + "source.p = source_mag * np.sin(2 * np.pi * source_freq * kgrid.t_array)\n", + "\n", + "# filter the source to remove high frequencies not supported by the grid\n", + "source.p = filter_time_series(kgrid, medium, source.p)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define simulation parameters and run simulations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# pre-allocate array for storing the output time series\n", + "single_element_data = np.zeros((Nt, points)) # noqa: F841\n", + "\n", + "# run a simulation for each of these sources to see the effect that the\n", + "# angle from the detector has on the measured signal\n", + "for source_loop in range(points):\n", + " # select a point source\n", + " source.p_mask = np.zeros(grid_size_points)\n", + " source.p_mask[unflatten_matlab_mask(source.p_mask, source_positions[source_loop] - 1)] = 1\n", + "\n", + " # run the simulation\n", + "\n", + " input_filename = f\"example_input_{source_loop + 1}_input.h5\"\n", + " pathname = gettempdir()\n", + " input_file_full_path = os.path.join(pathname, input_filename)\n", + " simulation_options = SimulationOptions(save_to_disk=True, input_filename=input_filename, data_path=pathname)\n", + " # run the simulation\n", + " sensor_data = kspaceFirstOrder2DC(\n", + " medium=medium,\n", + " kgrid=kgrid,\n", + " source=deepcopy(source),\n", + " sensor=deepcopy(sensor),\n", + " simulation_options=simulation_options,\n", + " execution_options=SimulationExecutionOptions(),\n", + " )\n", + " single_element_data[:, source_loop] = sensor_data['p'].sum(axis=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize recorded data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.imshow(circle + sensor.mask, \n", + " extent=[\n", + " kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, \n", + " kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], vmin=-1, vmax=1, cmap=get_color_map())\n", + "plt.ylabel('x-position [mm]')\n", + "plt.xlabel('y-position [mm]')\n", + "plt.axis('image')\n", + "\n", + "_, t_scale, t_prefix, _ = scale_SI(kgrid.t_array[-1])\n", + "\n", + "# Plot the time series recorded for each of the sources\n", + "plt.figure()\n", + "# Loop through each source and plot individually\n", + "for i in range(single_element_data.shape[1]):\n", + " plt.plot((kgrid.t_array * t_scale).squeeze(), single_element_data[:, i], label=f'Source {i+1}')\n", + "\n", + "plt.xlabel(f'Time [{t_prefix}s]')\n", + "plt.ylabel('Pressure [au]')\n", + "plt.title('Time Series For Each Source Direction')\n", + "\n", + "\n", + "# Calculate angle between source and center of detector face\n", + "angles = []\n", + "for source_position in source_positions:\n", + " x, y = ind2sub(kgrid.y.shape, source_position)\n", + " angles.append(np.arctan(kgrid.y[x, y] / kgrid.x[x, y]))\n", + "\n", + "# Plot the maximum amplitudes for each of the sources\n", + "plt.figure()\n", + "plt.plot(angles, np.max(single_element_data[200:350, :], axis=0), 'o', mfc='none')\n", + "plt.xlabel('Angle Between Source and Centre of Detector Face [rad]')\n", + "plt.ylabel('Maximum Detected Pressure [au]')\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "env_kwave", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/sd_directivity_modelling_2D/sd_directivity_modelling_2D.py b/examples/sd_directivity_modelling_2D/sd_directivity_modelling_2D.py new file mode 100644 index 00000000..0eef6c26 --- /dev/null +++ b/examples/sd_directivity_modelling_2D/sd_directivity_modelling_2D.py @@ -0,0 +1,127 @@ +import os +from copy import deepcopy +from tempfile import gettempdir + +import numpy as np +import matplotlib.pyplot as plt + +from kwave.data import Vector +from kwave.kgrid import kWaveGrid +from kwave.kmedium import kWaveMedium +from kwave.ksensor import kSensor +from kwave.ksource import kSource +from kwave.kspaceFirstOrder2D import kspaceFirstOrder2DC +from kwave.options.simulation_execution_options import SimulationExecutionOptions +from kwave.options.simulation_options import SimulationOptions +from kwave.utils.conversion import cart2grid +from kwave.utils.filters import filter_time_series +from kwave.utils.mapgen import make_cart_circle +from kwave.utils.matlab import ind2sub, matlab_find, unflatten_matlab_mask +from kwave.utils.colormap import get_color_map +from kwave.utils.data import scale_SI + + +grid_size_points = Vector([128, 128]) # [grid points] +grid_size_meters = Vector([50e-3, 50e-3]) # [m] +grid_spacing_meters = grid_size_meters / grid_size_points # [m] +kgrid = kWaveGrid(grid_size_points, grid_spacing_meters) + +# define the properties of the propagation medium +medium = kWaveMedium(sound_speed=1500) + +# define the array of time points [s] +Nt = 350 +dt = 7e-8 # [s] +kgrid.setTime(Nt, dt) + +sz = 20 # [grid points] +sensor_mask = np.zeros(grid_size_points) +sensor_mask[grid_size_points.x // 2, (grid_size_points.y // 2 - sz // 2) : (grid_size_points.y // 2 + sz // 2) + 1] = 1 +sensor = kSensor(sensor_mask) + + +# define equally spaced point sources lying on a circle centred at the +# centre of the detector face +radius = 30 # [grid points] +points = 11 +circle = make_cart_circle(radius * grid_spacing_meters.x, points, Vector([0, 0]), np.pi) + +# find the binary sensor mask most closely corresponding to the Cartesian +# coordinates from makeCartCircle +circle, _, _ = cart2grid(kgrid, circle) + +# find the indices of the sources in the binary source mask +source_positions = matlab_find(circle, val=1, mode="eq") + +source = kSource() +source_freq = 0.25e6 # [Hz] +source_mag = 1 # [Pa] +source.p = source_mag * np.sin(2 * np.pi * source_freq * kgrid.t_array) + +# filter the source to remove high frequencies not supported by the grid +source.p = filter_time_series(kgrid, medium, source.p) + +# pre-allocate array for storing the output time series +single_element_data = np.zeros((Nt, points)) # noqa: F841 + +# run a simulation for each of these sources to see the effect that the +# angle from the detector has on the measured signal +for source_loop in range(points): + # select a point source + source.p_mask = np.zeros(grid_size_points) + source.p_mask[unflatten_matlab_mask(source.p_mask, source_positions[source_loop] - 1)] = 1 + + # run the simulation + + input_filename = f"example_input_{source_loop + 1}_input.h5" + pathname = gettempdir() + input_file_full_path = os.path.join(pathname, input_filename) + simulation_options = SimulationOptions(save_to_disk=True, input_filename=input_filename, data_path=pathname) + # run the simulation + sensor_data = kspaceFirstOrder2DC( + medium=medium, + kgrid=kgrid, + source=deepcopy(source), + sensor=deepcopy(sensor), + simulation_options=simulation_options, + execution_options=SimulationExecutionOptions(), + ) + single_element_data[:, source_loop] = sensor_data["p"].sum(axis=1) + +plt.figure() +plt.imshow( + circle + sensor.mask, + extent=[kgrid.y_vec.min() * 1e3, kgrid.y_vec.max() * 1e3, kgrid.x_vec.max() * 1e3, kgrid.x_vec.min() * 1e3], + vmin=-1, + vmax=1, + cmap=get_color_map(), +) +plt.ylabel("x-position [mm]") +plt.xlabel("y-position [mm]") +plt.axis("image") + +_, t_scale, t_prefix, _ = scale_SI(kgrid.t_array[-1]) + +# Plot the time series recorded for each of the sources +plt.figure() +# Loop through each source and plot individually +for i in range(single_element_data.shape[1]): + plt.plot((kgrid.t_array * t_scale).squeeze(), single_element_data[:, i], label=f"Source {i+1}") + +plt.xlabel(f"Time [{t_prefix}s]") +plt.ylabel("Pressure [au]") +plt.title("Time Series For Each Source Direction") + + +# Calculate angle between source and center of detector face +angles = [] +for source_position in source_positions: + x, y = ind2sub(kgrid.y.shape, source_position) + angles.append(np.arctan(kgrid.y[x, y] / kgrid.x[x, y])) + +# Plot the maximum amplitudes for each of the sources +plt.figure() +plt.plot(angles, np.max(single_element_data[200:350, :], axis=0), "o", mfc="none") +plt.xlabel("Angle Between Source and Centre of Detector Face [rad]") +plt.ylabel("Maximum Detected Pressure [au]") +plt.show() diff --git a/kwave/__init__.py b/kwave/__init__.py index 43399a1f..4181acd5 100644 --- a/kwave/__init__.py +++ b/kwave/__init__.py @@ -7,9 +7,10 @@ import hashlib import json + # Test installation with: # python3 -m pip install -i https://test.pypi.org/simple/ --extra-index-url=https://pypi.org/simple/ k-Wave-python==0.3.0 -VERSION = "0.3.4" +VERSION = "0.3.5" # Constants and Configurations URL_BASE = "https://github.com/waltsims/" @@ -21,7 +22,9 @@ raise NotImplementedError(f"k-wave-python is currently unsupported on this operating system: {PLATFORM}.") # TODO: install directly in to /bin/ directory system directory is no longer needed +# TODO: depricate in 0.3.7 BINARY_PATH = Path(__file__).parent / "bin" / PLATFORM +BINARY_DIR = BINARY_PATH # add alias for BINARY_PATH for now WINDOWS_DLLS = [ diff --git a/kwave/executor.py b/kwave/executor.py index 8456fe5b..4834d8c9 100644 --- a/kwave/executor.py +++ b/kwave/executor.py @@ -6,11 +6,12 @@ import h5py import kwave +from kwave.options.simulation_execution_options import SimulationExecutionOptions from kwave.utils.dotdictionary import dotdict class Executor: - def __init__(self, execution_options, simulation_options): + def __init__(self, execution_options: SimulationExecutionOptions, simulation_options): self.execution_options = execution_options self.simulation_options = simulation_options @@ -21,15 +22,15 @@ def __init__(self, execution_options, simulation_options): self._make_binary_executable() def _make_binary_executable(self): - try: - self.execution_options.binary_path.chmod(self.execution_options.binary_path.stat().st_mode | stat.S_IEXEC) - except FileNotFoundError as e: + binary_path = self.execution_options.binary_path + if not binary_path.exists(): if kwave.PLATFORM == "darwin" and self.execution_options.is_gpu_simulation: raise ValueError( - "GPU simulations are currently not supported on MacOS. Try running the simulation on CPU by setting is_gpu_simulation=False." - ) from e - else: - raise e + "GPU simulations are currently not supported on MacOS. " + "Try running the simulation on CPU by setting is_gpu_simulation=False." + ) + raise FileNotFoundError(f"Binary not found at {binary_path}") + binary_path.chmod(binary_path.stat().st_mode | stat.S_IEXEC) def run_simulation(self, input_filename: str, output_filename: str, options: str): command = [str(self.execution_options.binary_path), "-i", input_filename, "-o", output_filename, options] diff --git a/kwave/options/simulation_execution_options.py b/kwave/options/simulation_execution_options.py index 3dcf42c6..7e46bc3c 100644 --- a/kwave/options/simulation_execution_options.py +++ b/kwave/options/simulation_execution_options.py @@ -1,120 +1,201 @@ -from dataclasses import dataclass -import os +from pathlib import Path from typing import Optional, Union +import os -from kwave import PLATFORM, BINARY_PATH +from kwave import PLATFORM, BINARY_DIR from kwave.ksensor import kSensor -@dataclass class SimulationExecutionOptions: - # Are we going to run the simulation on the GPU? - # Affects binary name and the way the simulation is run - is_gpu_simulation: bool = False - - # user defined location for the binary - binary_path: Optional[str] = BINARY_PATH - # user defined location for the binary - binary_name: Optional[str] = None - # user defined number of threads - kwave_function_name: Optional[str] = "kspaceFirstOrder3D" - # Whether to delete the input and output files after the simulation - delete_data: bool = True - # GPU device flag - device_num: Optional[int] = None - # number of threads - num_threads: Union[int, str] = "all" - - # user defined thread binding option - thread_binding: Optional[bool] = None - - # user input for system string - system_call: Optional[str] = None - verbose_level: int = 0 - - # determine whether chunking is handled automatically (the default), or manually - auto_chunking: Optional[bool] = True - - # show simulation log - show_sim_log: bool = True - - def __post_init__(self): - self.validate() - - if self.binary_name is None: + """ + A class to manage and configure the execution options for k-Wave simulations. + """ + + def __init__( + self, + is_gpu_simulation: bool = False, + binary_path: Optional[str] = None, + binary_dir: Optional[str] = None, + binary_name: Optional[str] = None, + kwave_function_name: Optional[str] = "kspaceFirstOrder3D", + delete_data: bool = True, + device_num: Optional[int] = None, + num_threads: Union[int, str] = "all", + thread_binding: Optional[bool] = None, + system_call: Optional[str] = None, + verbose_level: int = 0, + auto_chunking: Optional[bool] = True, + show_sim_log: bool = True, + ): + self.is_gpu_simulation = is_gpu_simulation + self._binary_path = binary_path + self._binary_name = binary_name + self._binary_dir = binary_dir + self.kwave_function_name = kwave_function_name + self.delete_data = delete_data + self.device_num = device_num + self.num_threads = num_threads + self.thread_binding = thread_binding + self.system_call = system_call + self.verbose_level = verbose_level + self.auto_chunking = auto_chunking + self.show_sim_log = show_sim_log + + @property + def num_threads(self) -> Union[int, str]: + return self._num_threads + + @num_threads.setter + def num_threads(self, value: Union[int, str]): + cpu_count = os.cpu_count() + if cpu_count is None: + raise RuntimeError("Unable to determine the number of CPUs on this system. Please specify the number of threads explicitly.") + + if value == "all": + value = cpu_count + + if not isinstance(value, int): + raise ValueError("Got {value}. Number of threads must be 'all' or a positive integer") + + if value <= 0 or value > cpu_count: + raise ValueError(f"Number of threads {value} must be a positive integer and less than total threads on the system {cpu_count}.") + + self._num_threads = value + + @property + def verbose_level(self) -> int: + return self._verbose_level + + @verbose_level.setter + def verbose_level(self, value: int): + if not (isinstance(value, int) and 0 <= value <= 2): + raise ValueError("Verbose level must be between 0 and 2") + self._verbose_level = value + + @property + def is_gpu_simulation(self) -> Optional[bool]: + return self._is_gpu_simulation + + @is_gpu_simulation.setter + def is_gpu_simulation(self, value: Optional[bool]): + "Set the flag to enable default GPU simulation. This option will supercede custom binary paths." + self._is_gpu_simulation = value + # Automatically update the binary name based on the GPU simulation flag + if value is not None: + self._binary_name = None + + @property + def binary_name(self) -> str: + valid_binary_names = ["kspaceFirstOrder-CUDA", "kspaceFirstOrder-OMP"] + if PLATFORM == "windows": + valid_binary_names = [name + ".exe" for name in valid_binary_names] + + if self._binary_name is None: + # set default binary name based on GPU simulation value + if self.is_gpu_simulation is None: + raise ValueError("`is_gpu_simulation` must be set to either True or False before determining the binary name.") + if self.is_gpu_simulation: - self.binary_name = "kspaceFirstOrder-CUDA" + self._binary_name = "kspaceFirstOrder-CUDA" else: - self.binary_name = "kspaceFirstOrder-OMP" + self._binary_name = "kspaceFirstOrder-OMP" - if PLATFORM == "windows": - if not self.binary_name.endswith(".exe"): - self.binary_name += ".exe" + if PLATFORM == "windows": + self._binary_name += ".exe" + elif self._binary_name not in valid_binary_names: + import warnings - self.binary_path = self.binary_path / self.binary_name + warnings.warn("Custom binary name set. Ignoring `is_gpu_simulation` state.") + return self._binary_name - def validate(self): - if isinstance(self.num_threads, int): - assert self.num_threads > 0 and self.num_threads != float("inf") - else: - assert self.num_threads == "all" - self.num_threads = None + @binary_name.setter + def binary_name(self, value: str): + self._binary_name = value - assert isinstance(self.verbose_level, int) and 0 <= self.verbose_level <= 2 + @property + def binary_path(self) -> Path: + if self._binary_path is not None: + return self._binary_path - def get_options_string(self, sensor: kSensor) -> str: - options_dict = {} - if self.device_num: - options_dict["-g"] = self.device_num + binary_dir = BINARY_DIR if self._binary_dir is None else self._binary_dir - if self.num_threads: - options_dict["-t"] = self.num_threads + if binary_dir is None: + raise ValueError("Binary directory is not specified.") - if self.verbose_level: - options_dict["--verbose"] = self.verbose_level + path = Path(binary_dir) / self.binary_name + if PLATFORM == "windows" and not path.name.endswith(".exe"): + path = path.with_suffix(".exe") + return path - options_string = "" - for flag, value in options_dict.items(): - options_string += f" {flag} {str(value)}" + @binary_path.setter + def binary_path(self, value: str): + # check if the binary path is a valid path + if not os.path.exists(value): + raise FileNotFoundError( + f"Binary path {value} does not exist. If you are trying to set `binary_dir`, use the `binary_dir` attribute instead." + ) + self._binary_path = value + + @property + def binary_dir(self) -> str: + return BINARY_DIR if self._binary_dir is None else self._binary_dir + + @binary_dir.setter + def binary_dir(self, value: str): + # check if binary_dir is a directory + if not os.path.isdir(value): + raise NotADirectoryError( + f"{value} is not a directory. If you are trying to set the `binary_path`, use the `binary_path` attribute instead." + ) + self._binary_dir = Path(value) + + def get_options_string(self, sensor: kSensor) -> str: + options_list = [] + if self.device_num is not None and self.device_num < 0: + raise ValueError("Device number must be non-negative") + if self.device_num is not None: + options_list.append(f" -g {self.device_num}") + + if self.num_threads is not None: + options_list.append(f" -t {self.num_threads}") + + if self.verbose_level > 0: + options_list.append(f" --verbose {self.verbose_level}") + + record_options_map = { + "p": "p_raw", + "p_max": "p_max", + "p_min": "p_min", + "p_rms": "p_rms", + "p_max_all": "p_max_all", + "p_min_all": "p_min_all", + "p_final": "p_final", + "u": "u_raw", + "u_max": "u_max", + "u_min": "u_min", + "u_rms": "u_rms", + "u_max_all": "u_max_all", + "u_min_all": "u_min_all", + "u_final": "u_final", + } - # check if sensor.record is given if sensor.record is not None: - record = sensor.record - - # set the options string to record the required output fields - record_options_map = { - "p": "p_raw", - "p_max": "p_max", - "p_min": "p_min", - "p_rms": "p_rms", - "p_max_all": "p_max_all", - "p_min_all": "p_min_all", - "p_final": "p_final", - "u": "u_raw", - "u_max": "u_max", - "u_min": "u_min", - "u_rms": "u_rms", - "u_max_all": "u_max_all", - "u_min_all": "u_min_all", - "u_final": "u_final", - } - for k, v in record_options_map.items(): - if k in record: - options_string = options_string + f" --{v}" - - if "u_non_staggered" in record or "I_avg" in record or "I" in record: - options_string = options_string + " --u_non_staggered_raw" - - if ("I_avg" in record or "I" in record) and ("p" not in record): - options_string = options_string + " --p_raw" + matching_keys = set(sensor.record).intersection(record_options_map.keys()) + for key in matching_keys: + options_list.append(f" --{record_options_map[key]}") + + if "u_non_staggered" in sensor.record or "I_avg" in sensor.record or "I" in sensor.record: + options_list.append(" --u_non_staggered_raw") + + if ("I_avg" in sensor.record or "I" in sensor.record) and ("p" not in sensor.record): + options_list.append(" --p_raw") else: - # if sensor.record is not given, record the raw time series of p - options_string = options_string + " --p_raw" + options_list.append(" --p_raw") - # check if sensor.record_start_imdex is given if sensor.record_start_index is not None: - options_string = options_string + " -s " + str(sensor.record_start_index) - return options_string + options_list.append(f" -s {sensor.record_start_index}") + + return " ".join(options_list) @property def env_vars(self) -> dict: diff --git a/pyproject.toml b/pyproject.toml index 6ca0d3c7..c9b8506c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,7 +8,7 @@ dynamic = ["version"] description = "Acoustics toolbox for time domain acoustic and ultrasound simulations in complex and tissue-realistic media." readme = "docs/README.md" license = { file = "LICENSE" } -requires-python = ">=3.9" +requires-python = ">=3.10" authors = [ { name = "Farid Yagubbayli", email = "farid.yagubbayli@tum.de" }, { name = "Walter Simson", email = "walter.simson@tum.de"} diff --git a/tests/test_executor.py b/tests/test_executor.py index f65f37d8..b3414fea 100644 --- a/tests/test_executor.py +++ b/tests/test_executor.py @@ -193,6 +193,7 @@ def test_executor_gpu_cuda_failure_darwin(self): # Configure the mock path object mock_binary_path = MagicMock(spec=Path) mock_binary_path.chmod.side_effect = FileNotFoundError + mock_binary_path.exists.return_value = False # Mock the execution options to use the mocked path mock_execution_options = MagicMock() diff --git a/tests/test_simulation_execution_options.py b/tests/test_simulation_execution_options.py index 2da09e72..efaeda0b 100644 --- a/tests/test_simulation_execution_options.py +++ b/tests/test_simulation_execution_options.py @@ -1,30 +1,85 @@ import unittest -from kwave.options.simulation_execution_options import SimulationExecutionOptions +from unittest.mock import Mock, patch from kwave.ksensor import kSensor -from unittest.mock import patch +from kwave import PLATFORM +from kwave.options.simulation_execution_options import SimulationExecutionOptions +import os + + +OMP_BINARY_NAME = "kspaceFirstOrder-OMP{}".format(".exe" if PLATFORM == "windows" else "") +CUDA_BINARY_NAME = "kspaceFirstOrder-CUDA{}".format(".exe" if PLATFORM == "windows" else "") class TestSimulationExecutionOptions(unittest.TestCase): def setUp(self): - # Set up a default kSensor object for testing - self.sensor = kSensor() - self.sensor.record = ["p", "u"] - self.sensor.record_start_index = 10 + self.default_options = SimulationExecutionOptions() + self.mock_sensor = Mock(spec=kSensor) + self.mock_sensor.record = ["p", "u_max"] + self.mock_sensor.record_start_index = 10 def test_default_initialization(self): - with patch("kwave.options.simulation_execution_options.PLATFORM", "linux"): - options = SimulationExecutionOptions() - self.assertFalse(options.is_gpu_simulation) - self.assertEqual(options.binary_name, "kspaceFirstOrder-OMP") - self.assertTrue(options.delete_data) - self.assertEqual(options.num_threads, None) # TODO: confusing logic here - self.assertEqual(options.verbose_level, 0) + """Test default values during initialization.""" + options = self.default_options + self.assertFalse(options.is_gpu_simulation) + self.assertIsNone(options._binary_path) + self.assertIsNone(options._binary_name) + self.assertEqual(options.kwave_function_name, "kspaceFirstOrder3D") + self.assertTrue(options.delete_data) + self.assertIsNone(options.device_num) + self.assertEqual(options.num_threads, os.cpu_count()) # "all" should default to CPU count + self.assertIsNone(options.thread_binding) + self.assertEqual(options.verbose_level, 0) + self.assertTrue(options.auto_chunking) + self.assertTrue(options.show_sim_log) - def test_gpu_simulation_initialization(self): - with patch("kwave.options.simulation_execution_options.PLATFORM", "linux"): - options = SimulationExecutionOptions(is_gpu_simulation=True) - self.assertTrue(options.is_gpu_simulation) - self.assertEqual(options.binary_name, "kspaceFirstOrder-CUDA") + def test_num_threads_setter_valid(self): + """Test setting a valid number of threads.""" + options = self.default_options + options.num_threads = os.cpu_count() + self.assertEqual(options.num_threads, os.cpu_count()) + + options.num_threads = "all" + self.assertEqual(options.num_threads, os.cpu_count()) + + def test_num_threads_setter_invalid(self): + """Test setting an invalid number of threads.""" + options = self.default_options + with self.assertRaises(ValueError): + options.num_threads = -1 + with self.assertRaises(ValueError): + options.num_threads = "invalid_value" + + def test_verbose_level_setter_valid(self): + """Test setting a valid verbose level.""" + options = self.default_options + for level in range(3): + options.verbose_level = level + self.assertEqual(options.verbose_level, level) + + def test_verbose_level_setter_invalid(self): + """Test setting an invalid verbose level.""" + options = self.default_options + with self.assertRaises(ValueError): + options.verbose_level = 3 + with self.assertRaises(ValueError): + options.verbose_level = -1 + + def test_is_gpu_simulation_setter(self): + """Test setting is_gpu_simulation and its impact on binary_name.""" + options = self.default_options + options.is_gpu_simulation = True + self.assertTrue(options.is_gpu_simulation) + self.assertEqual(options.binary_name, CUDA_BINARY_NAME) + + options.is_gpu_simulation = False + self.assertFalse(options.is_gpu_simulation) + self.assertEqual(options.binary_name, OMP_BINARY_NAME) + + def test_binary_name_custom(self): + """Test setting a custom binary name.""" + options = self.default_options + options.binary_name = "custom_binary" + self.assertEqual(options.binary_name, "custom_binary") def test_binary_name_extension_on_windows(self): with patch("kwave.options.simulation_execution_options.PLATFORM", "windows"): @@ -32,11 +87,25 @@ def test_binary_name_extension_on_windows(self): self.assertTrue(options.binary_name.endswith(".exe")) def test_get_options_string(self): - options = SimulationExecutionOptions() - options_string = options.get_options_string(self.sensor) - self.assertIn("--p_raw", options_string) - self.assertIn("--u_raw", options_string) - self.assertIn("-s 10", options_string) + """Test the get_options_string method with a mock sensor.""" + options = self.default_options + options.device_num = 1 + options.num_threads = os.cpu_count() + options.verbose_level = 2 + + options_string = options.get_options_string(self.mock_sensor) + expected_substrings = [" -g 1", f" -t {os.cpu_count()}", " --verbose 2", " --p_raw", " --u_max", " -s 10"] + for substring in expected_substrings: + self.assertIn(substring, options_string) + + def test_gpu_dependency_on_binary_name_and_path(self): + """Test that the binary_name and binary_path are updated correctly based on is_gpu_simulation.""" + options = SimulationExecutionOptions(is_gpu_simulation=True) + self.assertEqual(options.binary_name, CUDA_BINARY_NAME) + + options.is_gpu_simulation = False + self.assertEqual(options.binary_name, OMP_BINARY_NAME) + self.assertTrue(str(options.binary_path).endswith(OMP_BINARY_NAME)) def test_env_vars_linux(self): with patch("kwave.options.simulation_execution_options.PLATFORM", "linux"):