Optimisation toolbox

The toolbox connects a user’s gprMax model to an external optimiser. You define what may change and how to score a completed simulation. The toolbox proposes values, builds and runs the model, reads your score, and returns it to the optimiser. Geometry can change on every trial; fixed-geometry Studies remain a separate workflow.

The user guide lives beside the code in gprMax/toolboxes/Optimisation and is included in the main gprMax documentation. Start here; the campaign and developer guides are optional references for more advanced uses.

Start with a small working example

Use the Python environment in which this gprMax development checkout runs. From the checkout directory, install the starter’s optional optimisation package:

python -m pip install optuna==4.9.0
mkdir -p my_optimisation
cp gprMax/toolboxes/Optimisation/examples/start_here.py my_optimisation/my_model.py
python my_optimisation/my_model.py

The example creates a small dielectric-block model, measures the peak receiver field and adjusts permittivity to approach a target peak. Its twelve evaluations are a short demonstration, not a convergence guarantee. It saves results beside the copied script. Set RESULTS to a fresh directory before running again.

Copy the starter. Save callbacks in a Python file; workers must be able to import them. Keep the launch under if __name__ == "__main__": as shown in the example. Additional input files and imported helpers should be listed in files= so changes are tracked. Finish editing before launching a campaign.

The four sections you edit

Table 25 Your model file

Section	What you supply	Why it is needed
Parameters and targets	PARAMETERS and your target constants	Defines allowed trial values and the desired physical response.
Model	build_model(parameters)	Uses those values to return a complete, unexecuted gprMax.Scene.
Criterion	evaluate(parameters, output)	Reads and processes the result; returns one finite number to minimise.
Launch	optimise(...)	Chooses the optimiser, budget, resources and results directory.

The keys in PARAMETERS are your names. For example, parameters["length_mm"] is the proposed physical length, not an internal optimiser coordinate. The same dictionary reaches both user functions.

Parameter values and units

from gprMax.toolboxes.Optimisation import Real, Integer, Categorical

PARAMETERS = {
    "length_mm": Integer(26, 36, "mm", step=2),  # 26, 28, 30, 32, 34, 36
    "permittivity": Real(2.0, 8.0),
    "material_name": Categorical(("dielectric_a", "dielectric_b")),
}

Real uses continuous bounds; Integer uses inclusive integer bounds and an optional positive integer step; Categorical chooses one of the supplied values. Both integer endpoints must lie on the declared lattice. unit is a label: convert millimetres to metres yourself in the model with value * 1e-3. The toolbox does not change the mesh or interpret your parameter names.

Use TPE or RF for categorical parameters. GA, PSO and DE currently accept numeric parameters only. All five support integer steps. Do not round to a different geometry in your model without recording what was actually built. Choose a mesh and parameter increments that represent your intended dimensions.

Build and check one model first

Create the domain, mesh, materials, geometry, source, receivers and any required monitors in build_model. Return the Scene without calling gprMax.run. The toolbox runs it in an isolated worker. To check one concrete design, put this code in a saved script:

from gprMax.toolboxes.Optimisation import simulate
from gprMax.toolboxes.Optimisation.examples.start_here import build_model, evaluate

if __name__ == "__main__":
    values = {"permittivity": 4.0}
    output = simulate(model=build_model, parameters=values,
                      directory="results/one_model")
    print(output.file)          # Path to the solver's main HDF5 file.
    print(output.datasets())    # Dataset names that were actually saved.
    print(evaluate(values, output))

simulate takes concrete values, with no optimiser or parameter ranges. Verify the geometry and output before spending a large optimisation budget.

Choose the criterion and stopping condition

The line returned by evaluate defines the criterion. The optimiser always minimises: -30 is better than -20. To maximise a positive quantity such as gain, return its negative. To match a waveform, return a scalar error between simulated and measured traces. You may use NumPy, an existing postprocessor or your own functions; the optimiser never needs to understand the raw output.

These are separate settings:

• objective=evaluate supplies the function that computes the score.

• evaluations=120 limits candidate evaluations. Repeated candidates count. Population methods use whole batches and may leave a remainder unused.

• target_value=1.0 optionally stops when the best score is at most 1.0. All candidates in the current batch finish first. Omit it for a budget-only run. stop_on_target=False records target attainment while continuing the budget.

The target has the same units and sign as the returned score. For an S11 objective in dB, target_value=-30.0 means an acceptable S11 of -30 dB or better. Neither a used-up budget nor a flat best-score curve proves global optimality. Missing or invalid outputs raise an error; correct the model or evaluation before treating such cases as a deliberate physical penalty.

Run, inspect and change the optimiser

from gprMax.toolboxes.Optimisation import optimise
from gprMax.toolboxes.Optimisation.examples.start_here import PARAMETERS, build_model, evaluate

if __name__ == "__main__":
    result = optimise(
        parameters=PARAMETERS, model=build_model, objective=evaluate,
        directory="results/my_de_run", optimiser="de",
        evaluations=120, population_size=10, seed=7,
    )
    print(result.best_parameters)
    print(result.best_value)
    print(result.stop_reason)

Change optimiser to "tpe", "rf", "ga", "pso" or "de" after installing its package. population_size controls GA/PSO/DE; batch_size controls simultaneous proposals for TPE/RF. They are independent of the number of available worker processes.

Results and resources

directory is a fresh campaign directory, not an input-output folder you must populate yourself. The toolbox creates one candidate directory per evaluation, keeps its parameter record and solver log, and passes the correct result to your objective. output.file identifies that candidate’s main HDF5 output. result.trials records the evaluated parameters, scores and statuses; optimiser/result.json saves the search result. Processing artifacts saved with output.save_npz belong to the evaluation that produced them.

The default runs serially on the CPU. A local pool can run separate models in parallel, including one per GPU supported by your gprMax installation:

from gprMax.toolboxes.Optimisation import LocalPool

execution = LocalPool(workers=4, cpu_threads_per_worker=2)
# GPU alternative, with these device IDs visible to this process:
# execution = LocalPool(solver="cuda", devices=[0, 1], cpu_threads_per_worker=2)
# Add execution=execution to the optimise() call above.

Every worker builds a complete model; each needs enough memory. This distributes candidate simulations, not one FDTD domain. Use enough candidates per batch to fill the available workers. MPI, resource profiles, caching and native optimiser checkpoint recovery are covered in the execution guide. The simple optimise call creates a new campaign; use the explicit Campaign interface for checkpoint recovery and shared simulation caches.

External optimisation packages

Install packages in the same Python environment as gprMax. The versions below were used for the toolbox trials; they are tested versions, not a claim that they are the latest available releases. Optional packages are imported when their adapter is used.

Table 26 Tried integrations

Selection	Package tested	Public adapter	Method and parameter support
"tpe"	Optuna 4.9.0	OptunaTPE	Tree-structured Parzen Estimator; real, integer and categorical.
"rf"	scikit-optimize 0.10.2	SkoptRF	Random-forest surrogate with expected improvement by default; real, integer and categorical.
"ga"	pymoo 0.6.2	PymooGA	Genetic algorithm; real and integer.
"pso"	pymoo 0.6.2	PymooPSO	Particle swarm optimisation; real and integer.
"de"	pymoo 0.6.2	PymooDE	Differential evolution; real and integer.
Direct helper	SciPy 1.18.0	minimise_integer	A separate bounded scalar search for one integer variable.

For one package, use the appropriate command:

python -m pip install optuna==4.9.0
python -m pip install scikit-optimize==0.10.2
python -m pip install pymoo==0.6.2

To install the five general adapters together, use python -m pip install -r gprMax/toolboxes/Optimisation/requirements-optimisers.txt. SciPy is already a gprMax dependency; the direct helper is described separately below. Check the active interpreter with python -m pip show optuna pymoo scikit-optimize scipy if an import fails.

Switching packages without changing the problem

optimiser_choices.py imports the starter’s parameter definitions, builder and objective. Only the launch settings change. For example, from the checkout:

python -m gprMax.toolboxes.Optimisation.examples.optimiser_choices results/tpe --optimiser tpe
python -m gprMax.toolboxes.Optimisation.examples.optimiser_choices results/rf --optimiser rf
python -m gprMax.toolboxes.Optimisation.examples.optimiser_choices results/ga --optimiser ga
python -m gprMax.toolboxes.Optimisation.examples.optimiser_choices results/pso --optimiser pso
python -m gprMax.toolboxes.Optimisation.examples.optimiser_choices results/de --optimiser de

Run only the methods you want to try. Each command creates a separate twelve-trial search on the same small model. Population methods use six candidates per batch. These short runs demonstrate the connection to the solver; they do not rank the algorithms or establish convergence.

RF and TPE learn from evaluated parameter/score pairs. RF predicts promising candidates with its forest, but the objective sent back to the optimiser is still computed from the actual gprMax result. This is distinct from training a surrogate to replace the solver permanently.

GA, PSO and DE propose populations. Integer coordinates are mapped to the declared physical step before reaching the model. PSO retains its continuous internal positions and velocities; GA and DE repair their coordinates onto the allowed lattice. This implementation detail does not change the dictionary your model receives. Use the same parameter definitions and objective when comparing methods.

Setting algorithm options

Pass an adapter object instead of a string when you want package-specific settings. Use one of these objects as optimiser=adapter in the ordinary optimise call. Its seed and batch/population size take precedence over the string-selection settings:

from gprMax.toolboxes.Optimisation import OptunaTPE, SkoptRF, PymooGA, PymooPSO, PymooDE

adapter = OptunaTPE(seed=7, batch_size=2)
adapter = SkoptRF(seed=7, n_initial_points=8, batch_size=2)
adapter = PymooGA(seed=7, population_size=10)
adapter = PymooPSO(seed=7, population_size=10, restart_after=20)
adapter = PymooDE(seed=7, population_size=10, options={"variant": "DE/rand/1/bin"})

These are alternatives: choose one assignment. restart_after is an optional toolbox PSO policy that starts a fresh swarm after that many populations without improvement, while retaining the campaign’s best result and total budget. It is disabled by default. Native pymoo PSO defaults remain active unless overridden through options. A larger population or more evaluations can change results; a seed makes a particular configuration repeatable, not universally successful.

The official package documentation explains their algorithms and native options:

• Optuna TPE.

• scikit-optimize ask/tell Optimizer.

• pymoo: GA, PSO, DE.

SciPy: a separate one-variable example

scipy_integer.py uses minimise_integer to call the same small model at integer permittivities:

python -m gprMax.toolboxes.Optimisation.examples.scipy_integer results/scipy

This helper wraps SciPy’s bounded minimize_scalar, remembers scores for repeated integers, and checks endpoints and neighbours. It is intended for an approximately unimodal interval. Its callback owns the simulation calls and storage. maxiter limits SciPy’s internal iterations; explicit endpoint/neighbour checks are additional evaluations.

There is no optimiser="scipy" alias in the general framework. Use the five adapters above for several parameters, parallel populations, shared campaign caching or native optimiser checkpoints.

What has been tried

TPE, RF, GA, PSO and DE have been exercised on the thin-wire dipole, as well as on small integration problems. The direct SciPy integer helper was also tried on the dipole. PSO and DE both found the published Antenna I dimensions in our rectangular-patch adaptation. A short budget did not make every dipole method converge equally; the examples demonstrate reusable mechanics, not an optimiser leaderboard. See the patch comparison and its modelling limits below.

Reading and processing solver outputs

The objective receives evaluate(parameters, output) after a candidate finishes. output is a SimulationOutput object for that candidate; you do not search a result directory to guess which file belongs to it.

Table 27 Public access points

Interface	Result
output.file	Path to the main solver HDF5 file; pass it to h5py or another postprocessor.
output.datasets()	Dataset paths actually saved in that file; lists names without loading arrays.
output.receiver(name, component, grid="/")	One receiver trace, its sampling times and units.
output.port(name, grid="/")	Terminal S11 and input impedance spectra, frequency and validity masks.
output.save_npz(name, **arrays)	Save intermediate arrays for this evaluation, for example a filtered waveform.

The convenience readers live in quantities.py. results.py supplies the methods above; processing.py manages evaluation records and saved artifacts. The same readers work without a campaign:

from gprMax.toolboxes.Optimisation import read_receiver, read_port

trace = read_receiver("output.h5", "probe", "Ez")
spectrum = read_port("output.h5", "frills/frill1")

These two calls illustrate different output types: your model must actually have the named receiver or terminal. The directory name does not determine which physical quantities exist.

ReceiverTrace fields

output.receiver("probe", "Ez") selects the receiver ID from your model and one saved component. It returns materialised arrays; the reader closes the file before returning.

Table 28 ReceiverTrace

Field	Meaning
values	One-dimensional sample array for the selected component.
time	Matching one-dimensional time array, in seconds.
unit	Component unit: for example V/m for E, A/m for H, or A for current.
dataset	Exact HDF5 dataset path used by the reader.
source	Source HDF5 file path.

Sampling requires the dataset’s explicit sample interval and offset metadata. Electric and magnetic quantities need not have identical sample times. Compare trace.time arrays before subtracting signals; the waveform example shows alignment and unit checks. This reader handles one-dimensional traces, not a merged B-scan matrix. Use grid to select an explicit subgrid rather than assuming its sampling matches the main grid.

PortSpectrum fields

Terminal ports, transmission lines and magnetic frills are selected under ports, tls and frills respectively. Prefer a qualified name such as "ports/feed" or "frills/frill1". An unqualified ID works only if unique. The name must match the saved solver output; the first magnetic frill is normally frill1. Inspect the file if unsure.

Table 29 PortSpectrum

Field	Meaning
frequency	One-dimensional frequency samples in Hz.
s11	Complex reflection coefficient, dimensionless; not dB.
valid	Boolean S11 validity mask with the same shape as frequency.
impedance	Complex input impedance in ohms.
impedance_valid	Separate Boolean impedance mask; do not substitute the S11 mask.
reference_impedance	Reference impedance in ohms used for the saved terminal spectrum.
tail_relative_db	Saved time-tail diagnostic in dB; NaN if absent.
independent_frequency_resolution_hz	Saved independent frequency resolution, or None if unavailable.
source / group	File path and exact terminal group path.

spectrum.at(frequency_hz) returns complex S11 using adjacent valid frequency samples. Interpolation acts on its real and imaginary parts. It refuses invalid gaps and extrapolation; it does not improve independent frequency resolution. Convert to dB explicitly in your criterion:

import numpy as np

def evaluate(parameters, output):
    spectrum = output.port("frills/frill1")
    reflection = spectrum.at(3.1e9)
    return float(20 * np.log10(max(abs(reflection), 1e-12)))

This minimises S11 in dB at one frequency. It does not locate the deepest dip anywhere in the spectrum. To plot or reduce a band, select valid samples inside your band first. Never replace invalid bins with zero: zero reflection would look like a perfect match. Missing groups, datasets or required attributes raise KeyError. A missing or ambiguous unqualified name, malformed data or invalid interpolation raises ValueError.

Voltages, currents and other outputs

The terminal HDF5 groups can contain V/I time traces and spectra in addition to S11 and Zin. PortSpectrum is a compact S11/Zin view, not a complete terminal record. For the existing broader terminal reader, use:

from gprMax.toolboxes.Plotting.plot_port import discover_port_outputs, read_port_output

print(discover_port_outputs(output.file))
port = read_port_output(output.file, "frills/frill1")
for trace in port.time_traces:
    print(trace.name, trace.quantity, trace.time.shape, trace.values.shape)

That reader lives in gprMax/toolboxes/Plotting/plot_port.py. Its PortData contains time_traces, spectral_traces, s_parameters, impedances, admittances, validity_masks, diagnostics and metadata. Depending on the terminal type and solver output, traces can include incident and total V/I, generator voltage and network or loop current. Use each time trace’s own time array; source-specific V and I traces can differ in length and time offset.

Table 30 Current coverage in the optimisation toolbox

Output family	How to use it now
Receivers	receiver / read_receiver for one saved component trace.
Source terminals	port / read_port for S11/Zin; the broader Plotting reader for V/I and other saved terminal fields.
Other source records	Inspect and read the actual HDF5 datasets through output.file.
Eigenmode ports	Read modal data with a source-specific postprocessor or h5py; there is no dedicated modal reader in Optimisation yet.
Antenna/scattering metrics	Configure the required monitors and postprocessing in your model/workflow, then reduce their results to your objective.

Gain, directivity, patterns, efficiency, RCS and SAR are not automatically computed by this toolbox. Their definition and required fields depend on the experiment. The objective may call an existing calculation using output.file or read other artifacts your workflow explicitly creates. Raw h5py access exposes stored data; interpreting axes, normalisation and units remains the postprocessor’s responsibility. A dedicated convenience reader for every output family has not yet been implemented.

Inspect first, then write the calculation

inspect_output.py prints group attributes and dataset names, shapes, types and attributes without loading all arrays. It can also demonstrate the public receiver/terminal reader:

python -m gprMax.toolboxes.Optimisation.examples.inspect_output path/to/output.h5
python -m gprMax.toolboxes.Optimisation.examples.inspect_output path/to/output.h5 --group frills
python -m gprMax.toolboxes.Optimisation.examples.inspect_output path/to/output.h5 --terminal frills/frill1 --frequency 3.1e9

Replace the path and selection with your actual file and source. For custom data, copy an observed dataset path rather than assuming names or units:

import h5py

with h5py.File(output.file, "r") as handle:
    print(list(handle.keys()))
    # values = handle["/the/dataset/path/you/inspected"][...]

Use output.save_npz("processed.npz", time_s=trace.time, signal=trace.values) to save your intermediate calculation in the objective. Names such as time_s and signal are yours, not mandatory toolbox fields. Optional units={"time_s": "s", "signal": "V/m"} records units with the processing artifact. Return the final finite scalar after any filtering, transformation or comparison. The optimiser interface stays the same.

Reusable examples and antenna validation

The main examples use the same four sections: parameters, model, objective and launch. Larger files put supporting checks after the main workflow and signpost each call. Choose the closest example, copy it, and edit the first three sections to define your own experiment.

Table 31 Which example to copy

File	Purpose
start_here.py	Small dielectric block; one parameter and a receiver-peak objective.
optimiser_choices.py	Change the optimiser while importing the same starter model and objective.
waveform_matching.py	Two parameters; load reference data, align waveforms and minimise normalised RMS error.
thin_wire_dipole.py	Change a symmetric dipole length so its S11 dip approaches 1 GHz.
rectangular_patch.py	Change patch length, width and feed offset; minimise S11 at 3.1 GHz.
inspect_output.py	Explore an existing output file and demonstrate public readers; no simulation.
scipy_integer.py	Direct SciPy helper for a single bounded integer.

dielectric_block.py, compare_dipole_optimisers.py and the advanced examples demonstrate explicit campaign/benchmark machinery. Start with the files above before using those additional interfaces.

Waveform matching: custom processing

Copy waveform_matching.py and reference_waveform.npz to the same folder. The bundled data are synthetic, generated with width 12 mm and relative permittivity 4. See reference notes.

Replace the reference file and loading block for measured data. Check time and amplitude units, align sampling times and define the error calculation in evaluate. The example saves comparison.npz for every scored trial so you can inspect the reference, simulated trace and residual. The optimiser receives only the normalised RMS error. Width is rounded to the example’s 2 mm mesh by its model builder; permittivity is continuous.

Thin-wire dipole: resonance location

thin_wire_dipole.py uses two wire arms and a one-cell voltage-source feed gap. The adjustable arm_cells is 54–83 cells per arm. With fixed DZ_M=0.001, total length is (2 * arm_cells + 1) * DZ_M: 109–167 mm in symmetric 2 mm increments.

The objective is the absolute difference, in Hz, between the estimated S11-dip frequency and 1 GHz. The helper at the end checks the spectrum and fits the dip between bins. This differs from minimising S11 directly at 1 GHz. The 100 ns time window gives approximately 10 MHz independent frequency resolution; interpolation is an estimate within that simulation, not a finer mesh or a longer recording. If you refine dz, review the cell-count bounds to preserve the physical lengths you intend to search.

Rectangular patch: literature-based validation

rectangular_patch.py is based on N. Jin and Y. Rahmat-Samii, “Parallel particle swarm optimization and finite-difference time-domain (PSO/FDTD) algorithm for multiband and wide-band patch antenna designs,” IEEE Transactions on Antennas and Propagation, 53(11), 3459–3468, 2005. Paper DOI. The example follows the single-frequency rectangular-patch case in Section III and compares against Antenna I in Table I. It is an adaptation with an equivalent magnetic-frill feed.

The fixed board/ground is 60 × 60 mm, substrate height 3 mm and relative permittivity 2.2. The mesh is 1 mm and the time window 100 ns. A 3 mm thin-wire probe spans the substrate, with a magnetic-frill source representing the coaxial aperture at the ground. Its assumed radius is 0.2 mm and reference impedance 50 ohms. These feed assumptions are documented in the source and patch notes.

Table 32 Physical parameter ranges, in millimetres

Parameter	Paper’s stated space	Reusable example
Patch length L	0 < L < 44	26, 28, 30, 32, 34, 36
Patch width W	0 < W < 44	12, 14, …, 44
Feed offset x	0 < x < 22; x < L/2	1, 2, …, 12; always inside the patch

The feed offset is measured from the patch centre, not the edge. Positive x in the parameter dictionary moves the feed toward the model’s negative x axis. Even L/W keep the centre and centreline on the fixed design lattice. These are physical millimetre increments; Integer(..., step=2) expresses them through the public interface for all five optimisers. This is a restricted benchmark space, not the paper’s full search.

The criterion is 20 * log10(abs(S11(3.1 GHz))). Lower is better. The paper’s single-frequency fitness adds 50 to this dB value; that constant does not change which candidate minimises it. The example does not optimise the broadband Antenna II criterion.

Both tested searches found L = 30 mm, W = 18 mm, x = 4 mm, the Table I Antenna I dimensions. The independently simulated reference and the selected result in this model give S11 = -31.5375 dB at 3.1 GHz; the paper reports approximately -28 dB. The agreement is in the design dimensions, not an exact reproduction of the paper’s feed, spectrum or PSO implementation.

Table 33 Recorded runs, seed 7 and ten candidates per population

Method	First reference design	Total proposals	New simulations	Stop
PSO	Population 31, proposal 304	550	68	User stopped after population 55
DE	Population 47, proposal 467	470	100	Reference score reached

The reference geometry was not seeded. Cached outputs were reused only after a candidate was proposed; counts of new simulations therefore differ from proposal counts. PSO used an optional restart after 20 stagnant populations, which first occurred after population 51, after finding the reference. DE used DE/rand/1/bin. The two single-seed runs validate the connection between the optimisers, geometry and objective; they do not establish an algorithm ranking.

Fig. 60 Recorded benchmark runs. The experiment’s cached runs and stopping settings differ from the short default example.

For a fresh DE run, save this as a script and use a new directory:

from gprMax.toolboxes.Optimisation import LocalPool
from gprMax.toolboxes.Optimisation.examples.rectangular_patch import run_example

if __name__ == "__main__":
    run_example(
        "results/patch_de", optimiser="de", evaluations=2000,
        population_size=10, seed=7, target_value=-30.0,
        execution=LocalPool(workers=2, cpu_threads_per_worker=2),
    )

This first checks one starting design, then starts the search. The starting simulation is additional to the optimisation budget and does not seed the optimiser. Change "de" to "pso" or another supported method. Set target_value=None to use the evaluation limit alone. The source file’s default budget is only 40 proposals; it is a short trial and need not find the reference. Fresh 100 ns simulations can make longer runs take hours.

The reusable interface does not preload the experiment’s private caches or resume its checkpoints. The model still needs mesh/time-window and feed checks for any new antenna application. This benchmark validates framework operation on this model; it is not a mesh-convergence study.

Further reference

• Campaign API: multiple scenarios, builders and stored records.

• Code walkthrough: follows a candidate from proposal to score.

• Developer guide: module responsibilities and extension points.

• Adapter contract: add another optimiser without changing model code.