Collect & Compare Example Inputs/Outputs

Author: farid-nucsai

Makefile

@@ -10,14 +10,29 @@ KWAVE_MATLAB_PATH = $(abspath ../k-wave) # Get absolute path of k-wave directory
 
 # Define the artifact directory
 COLLECTED_VALUES_DIR = $(abspath tests/matlab_test_data_collectors/python_testers/collectedValues)
+EXAMPLE_OUTPUT_DIR = /tmp/example_runs
+COMPARE_OUTPUT_DIR_A ?= /tmp/example_runs
+COMPARE_OUTPUT_DIR_B ?= /tmp/example_runs_1
 
 # Default target
 all: run-examples test
 
 # Target to run all examples
 run-examples:
 	@echo "Running all examples..."
-	@MPLBACKEND=Agg $(PYTHON) run_examples.py
+	@PYTHONPATH="$(CURDIR):$$PYTHONPATH" MPLBACKEND=Agg $(PYTHON) run_examples.py
+
+# Target to run all examples and collect non-dated HDF5 inputs/outputs
+run-examples-no-dates:
+	@echo "Running all examples with non-dated HDF5 filenames..."
+	@rm -rf $(EXAMPLE_OUTPUT_DIR)
+	@PYTHONPATH="$(CURDIR):$$PYTHONPATH" KWAVE_USE_DATED_FILENAMES=0 MPLBACKEND=Agg $(PYTHON) run_examples.py
+	@echo "Collected example files in $(EXAMPLE_OUTPUT_DIR)"
+
+# Target to compare two example output roots
+compare-example-outputs:
+	@echo "Comparing HDF outputs: $(COMPARE_OUTPUT_DIR_A) vs $(COMPARE_OUTPUT_DIR_B)"
+	@PYTHONPATH="$(CURDIR):$$PYTHONPATH" $(PYTHON) scripts/compare_example_outputs.py $(COMPARE_OUTPUT_DIR_A) $(COMPARE_OUTPUT_DIR_B)
 
 # Target to run pytest, which depends on running the MATLAB script first
 test: $(COLLECTED_VALUES_DIR)
@@ -41,4 +56,4 @@ clean-collected_values:
 	@echo "Cleaning collected values directory..."
 	@rm -rf $(COLLECTED_VALUES_DIR)
 
-.PHONY: all run-examples run-tests clean-python clean-collected_values clean
+.PHONY: all run-examples run-examples-no-dates compare-example-outputs run-tests clean-python clean-collected_values clean

examples/checkpointing/checkpoint.py

@@ -1,6 +1,5 @@
-from copy import copy
 from pathlib import Path
-from tempfile import TemporaryDirectory
+from shutil import copy2
 
 import numpy as np
 
@@ -15,6 +14,8 @@
 from kwave.utils.filters import smooth
 from kwave.utils.mapgen import make_ball
 
+EXAMPLE_OUTPUT_DIR = Path("/tmp/example_runs/checkpointing")
+
 # This script demonstrates how to use the checkpointing feature of k-Wave.
 # It runs the same simulation twice, with the second run starting from a
 # checkpoint file created during the first run.
@@ -30,6 +31,17 @@
 # Note: checkpoint timesteps and checkpoint interval must be integers.
 
 
+def _next_available_path(path: Path) -> Path:
+    if not path.exists():
+        return path
+    suffix = 1
+    while True:
+        candidate = path.with_name(f"{path.stem}_{suffix}{path.suffix}")
+        if not candidate.exists():
+            return candidate
+        suffix += 1
+
+
 def make_simulation_parameters(directory: Path, checkpoint_timesteps: int):
     """
     See the 3D FFT Reconstruction For A Planar Sensor example for context.
@@ -71,6 +83,7 @@ def make_simulation_parameters(directory: Path, checkpoint_timesteps: int):
         pml_inside=False,
         smooth_p0=False,
         data_cast="single",
+        allow_file_overwrite=True,
         input_filename=input_filename,
         output_filename=output_filename,
     )
@@ -84,23 +97,31 @@ def main():
     # the halfway point.
     checkpoint_timesteps = 106
 
-    with TemporaryDirectory() as tmpdir:
-        tmpdir = Path(tmpdir)
-        # create the simulation parameters for the 1st run
-        kgrid, medium, source, sensor, simulation_options, execution_options = make_simulation_parameters(
-            directory=tmpdir, checkpoint_timesteps=checkpoint_timesteps
-        )
-        kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options)
-        print("Temporary directory contains 3 files (including checkpoint):")
-        print("\t-", "\n\t- ".join([f.name for f in tmpdir.glob("*.h5")]))
-
-        # create the simulation parameters for the 2nd run
-        kgrid, medium, source, sensor, simulation_options, execution_options = make_simulation_parameters(
-            directory=tmpdir, checkpoint_timesteps=checkpoint_timesteps
-        )
-        kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options)
-        print("Temporary directory contains 2 files (checkpoint has been deleted):")
-        print("\t-", "\n\t- ".join([f.name for f in tmpdir.glob("*.h5")]))
+    EXAMPLE_OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
+    for filename in ("kwave_input.h5", "kwave_output.h5", "kwave_checkpoint.h5"):
+        (EXAMPLE_OUTPUT_DIR / filename).unlink(missing_ok=True)
+
+    # create the simulation parameters for the 1st run
+    kgrid, medium, source, sensor, simulation_options, execution_options = make_simulation_parameters(
+        directory=EXAMPLE_OUTPUT_DIR,
+        checkpoint_timesteps=checkpoint_timesteps,
+    )
+    kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options)
+    copy2(simulation_options.input_filename, _next_available_path(EXAMPLE_OUTPUT_DIR / "run_1_kwave_input.h5"))
+    copy2(simulation_options.output_filename, _next_available_path(EXAMPLE_OUTPUT_DIR / "run_1_kwave_output.h5"))
+    print("Checkpoint output directory after run_1:")
+    print("\t-", "\n\t- ".join([f.name for f in EXAMPLE_OUTPUT_DIR.glob("*.h5")]))
+
+    # create the simulation parameters for the 2nd run
+    kgrid, medium, source, sensor, simulation_options, execution_options = make_simulation_parameters(
+        directory=EXAMPLE_OUTPUT_DIR,
+        checkpoint_timesteps=checkpoint_timesteps,
+    )
+    kspaceFirstOrder3D(kgrid, source, sensor, medium, simulation_options, execution_options)
+    copy2(simulation_options.input_filename, _next_available_path(EXAMPLE_OUTPUT_DIR / "run_2_kwave_input.h5"))
+    copy2(simulation_options.output_filename, _next_available_path(EXAMPLE_OUTPUT_DIR / "run_2_kwave_output.h5"))
+    print("Checkpoint output directory after run_2:")
+    print("\t-", "\n\t- ".join([f.name for f in EXAMPLE_OUTPUT_DIR.glob("*.h5")]))
 
 
 if __name__ == "__main__":

examples/sd_directivity_modelling_2D/sd_directivity_modelling_2D.py

@@ -1,6 +1,4 @@
-import os
 from copy import deepcopy
-from tempfile import gettempdir
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -73,9 +71,8 @@
     # 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)
+    output_filename = f"example_output_{source_loop + 1}_output.h5"
+    simulation_options = SimulationOptions(save_to_disk=True, input_filename=input_filename, output_filename=output_filename)
     # run the simulation
     sensor_data = kspaceFirstOrder2DC(
         medium=medium,

examples/sd_focussed_detector_2D/sd_focussed_detector_2D.py

@@ -1,9 +1,7 @@
 # # Focussed Detector In 2D Example
 # This example shows how k-Wave-python can be used to model the output of a focused semicircular detector, where the directionality arises from spatially averaging across the detector surface. Unlike the original example in k-Wave, this example does not visualize the simulation, as this functionality is not intrinsically supported by the accelerated binaries.
 
-import os
 from copy import deepcopy
-from tempfile import gettempdir
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -43,9 +41,7 @@
 
 # ## Define simulation parameters
 input_filename = "example_sd_focused_2d_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)
+simulation_options = SimulationOptions(save_to_disk=True, input_filename=input_filename)
 
 
 # ## Run simulation with first source

examples/sd_focussed_detector_3D/sd_focussed_detector_3D.py

@@ -1,9 +1,7 @@
 # Focussed Detector In 3D Example
 # This example shows how k-Wave can be used to model the output of a focussed bowl detector where the directionality arises from spatially averaging across the detector surface.
 
-import os
 from copy import deepcopy
-from tempfile import gettempdir
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -33,10 +31,8 @@
 _ = kgrid.makeTime(medium.sound_speed)
 
 input_filename = "example_sd_focused_3d_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, pml_size=10, data_cast="single"
+    save_to_disk=True, input_filename=input_filename, pml_size=10, data_cast="single"
 )
 
 # create a concave sensor

examples/us_bmode_linear_transducer/us_bmode_linear_transducer.py

@@ -25,7 +25,7 @@
 
 # simulation settings
 DATA_CAST = "single"
-RUN_SIMULATION = False
+RUN_SIMULATION = True
 
 pml_size_points = Vector([20, 10, 10])  # [grid points]
 grid_size_points = Vector([256, 128, 128]) - 2 * pml_size_points  # [grid points]

examples/us_bmode_phased_array/us_bmode_phased_array.py

@@ -22,6 +22,7 @@
 # simulation settings
 DATA_CAST = "single"
 RUN_SIMULATION = True
+RNG_SEED = 123456
 
 
 pml_size_points = Vector([15, 10, 10])  # [grid points]
@@ -83,18 +84,19 @@
 
 not_transducer = NotATransducer(transducer, kgrid, **not_transducer)
 
+rng = np.random.default_rng(RNG_SEED)
 
 # Define a random distribution of scatterers for the medium
 background_map_mean = 1
 background_map_std = 0.008
-background_map = background_map_mean + background_map_std * np.random.randn(kgrid.Nx, kgrid.Ny, kgrid.Nz)
+background_map = background_map_mean + background_map_std * rng.standard_normal((kgrid.Nx, kgrid.Ny, kgrid.Nz))
 
 sound_speed_map = c0 * background_map
 density_map = rho0 * background_map
 
 
 # Define a random distribution of scatterers for the highly scattering region
-scattering_map = np.random.randn(kgrid.Nx, kgrid.Ny, kgrid.Nz)
+scattering_map = rng.standard_normal((kgrid.Nx, kgrid.Ny, kgrid.Nz))
 scattering_c0 = np.clip(c0 + 25 + 75 * scattering_map, 1400, 1600)
 scattering_rho0 = scattering_c0 / 1.5
 

kwave/kWaveSimulation_helper/save_to_disk_func.py

@@ -15,7 +15,14 @@
 
 
 def save_to_disk_func(
-    kgrid: kWaveGrid, medium: kWaveMedium, source, opt: SimulationOptions, auto_chunk: bool, values: dotdict, flags: dotdict
+    kgrid: kWaveGrid,
+    medium: kWaveMedium,
+    source,
+    opt: SimulationOptions,
+    auto_chunk: bool,
+    values: dotdict,
+    flags: dotdict,
+    input_filename: str | None = None,
 ):
     # update command line status
     logging.log(logging.INFO, f"  precomputation completed in {scale_time(TicToc.toc())}")
@@ -63,7 +70,7 @@ def save_to_disk_func(
     # =========================================================================
 
     remove_z_dimension(float_variables, kgrid.dim)
-    save_file(opt.input_filename, integer_variables, float_variables, opt.hdf_compression_level, auto_chunk=auto_chunk)
+    save_file(input_filename or opt.input_filename, integer_variables, float_variables, opt.hdf_compression_level, auto_chunk=auto_chunk)
 
     # update command line status
     logging.log(logging.INFO, f"  completed in {scale_time(TicToc.toc())}")

kwave/kspaceFirstOrder2D.py

@@ -12,7 +12,7 @@
 from kwave.kWaveSimulation import kWaveSimulation
 from kwave.kWaveSimulation_helper import retract_transducer_grid_size, save_to_disk_func
 from kwave.options.simulation_execution_options import SimulationExecutionOptions
-from kwave.options.simulation_options import SimulationOptions
+from kwave.options.simulation_options import SimulationOptions, resolve_filenames_for_run
 from kwave.utils.dotdictionary import dotdict
 from kwave.utils.interp import interpolate2d
 from kwave.utils.pml import get_pml
@@ -300,6 +300,7 @@ def kspaceFirstOrder2D(
     if options.save_to_disk:
         # store the pml size for resizing transducer object below
         retract_size = [[options.pml_x_size, options.pml_y_size, options.pml_z_size]]
+        input_filename, output_filename = resolve_filenames_for_run(k_sim.options)
 
         # run subscript to save files to disk
         save_to_disk_func(
@@ -348,6 +349,7 @@ def kspaceFirstOrder2D(
                     "cuboid_corners": k_sim.cuboid_corners,
                 }
             ),
+            input_filename=input_filename,
         )
 
         # run subscript to resize the transducer object if the grid has been expanded
@@ -359,5 +361,5 @@ def kspaceFirstOrder2D(
 
         executor = Executor(simulation_options=simulation_options, execution_options=execution_options)
         executor_options = execution_options.as_list(sensor=k_sim.sensor)
-        sensor_data = executor.run_simulation(k_sim.options.input_filename, k_sim.options.output_filename, options=executor_options)
+        sensor_data = executor.run_simulation(input_filename, output_filename, options=executor_options)
         return sensor_data

kwave/kspaceFirstOrder3D.py

@@ -12,7 +12,7 @@
 from kwave.kWaveSimulation import kWaveSimulation
 from kwave.kWaveSimulation_helper import retract_transducer_grid_size, save_to_disk_func
 from kwave.options.simulation_execution_options import SimulationExecutionOptions
-from kwave.options.simulation_options import SimulationOptions
+from kwave.options.simulation_options import SimulationOptions, resolve_filenames_for_run
 from kwave.utils.dotdictionary import dotdict
 from kwave.utils.interp import interpolate3d
 from kwave.utils.pml import get_pml
@@ -313,6 +313,7 @@ def kspaceFirstOrder3D(
     if options.save_to_disk:
         # store the pml size for resizing transducer object below
         retract_size = [[options.pml_x_size, options.pml_y_size, options.pml_z_size]]
+        input_filename, output_filename = resolve_filenames_for_run(k_sim.options)
 
         # run subscript to save files to disk
         save_to_disk_func(
@@ -361,6 +362,7 @@ def kspaceFirstOrder3D(
                     "cuboid_corners": k_sim.cuboid_corners,
                 }
             ),
+            input_filename=input_filename,
         )
 
         # run subscript to resize the transducer object if the grid has been expanded
@@ -372,5 +374,5 @@ def kspaceFirstOrder3D(
 
         executor = Executor(simulation_options=simulation_options, execution_options=execution_options)
         executor_options = execution_options.as_list(sensor=k_sim.sensor)
-        sensor_data = executor.run_simulation(k_sim.options.input_filename, k_sim.options.output_filename, options=executor_options)
+        sensor_data = executor.run_simulation(input_filename, output_filename, options=executor_options)
         return sensor_data