Sub-microsecond ADC delays are dropped without notification when exporting to .seq file

[jweine]

When adding ADC events that are starting on a valid time according to the adc-raster time but is smaller than a micro second the delay is simply dropped when exported to a .seq file.

This consistent with the file in section 2.8.3, as the table states delays are specified in integers of microseconds.

Below I have a minimally working example that creates a readout trapezoidal and a corresponding ADC from a number of samples and a flat-duration. To obtain a valid dwell time, the actual ADC duration is recomputed to be slighlty shorter, and subsequently the adc is shifted to be centered with respect to the flat-top duration.

Expected behaviour
I suspect severity of this depends on the choses raster_times, and in most cases with the default siemens raster times the differences might be neglegible. However, a Userwarning that states that this conversion is happening might be a good idea, so the user can choose e.g. a different number of samples or delibarately choose to ignore this difference.

To Reproduce

import pypulseq as pp
import numpy as np
system = pp.Opts(max_grad=32, grad_unit="mT/m", max_slew=130, slew_unit="T/m/s", 
                 rf_ringdown_time=30e-6, rf_dead_time=100e-6)

# Showcase invalid dwell time 
Nx = 121
flat_time = 1e-3
readout_time = flat_time
dwell_time = readout_time / Nx
print(f"Not a valid dwell time: {dwell_time}")
>>> Not a valid dwell time: 8.264462809917356e-06

# Define valid dwell time
Nx = 121
flat_time = 1e-3
dwell = np.floor(flat_time / Nx / system.adc_raster_time) * system.adc_raster_time
readout_time = dwell * Nx
difference = flat_time - readout_time
print(f"Difference ADC/Flat top: {difference}")
>>> Difference ADC/Flat top: 7.800000000000081e-06

# Center the adc on flat-top duration to ensure symmetric k-space sampling
adc_delay = difference / 2
print(f"Padding prior to ADC: {adc_delay}")
>>> Padding prior to ADC: 3.900000000000041e-06

fov = 220e-3
delta_k = 1 / fov
k_width = Nx * delta_k

seq = pp.Sequence()
gx = pp.make_trapezoid(channel="x", system=system,
                       amplitude=k_width / readout_time, flat_time=flat_time)
adc = pp.make_adc(num_samples=Nx, duration=readout_time, delay=adc_delay)
seq.add_block(gx, adc)
seq.write("dummy_adc.seq")

Which results in the following file, stating a 4us delay instead of a 3.9us delay.

[btasdelen]

Hi @jweine. As Siemens interpreter also uses microsecond for that delay, as you said this will be consistent with Siemens convention. In Siemens, as the time is an integer, that (flat_time - readout_time)/2 division will silently do the same thing.

We may be able to add a check and print a warning if the difference after integer cast is larger than eps.

[FrankZijlstra]

The functionality to check the timing of your sequence already exists:

    ok, error_report = seq.check_timing()
    if ok:
        print("Timing check passed successfully")
    else:
        print("Timing check failed. Error listing follows:")
        [print(e) for e in error_report]

This works for any timing related value, e.g. gradient timing as well. You should really always use this. It might make sense for seq.write to include a timing check.

I did notice there's a little bug in check_timing in that it does not print multiple timing errors within the same block (e.g. in the example if you provide a bad flat_time, it will still only report the ADC delay time error).

I've not been happy about the way to print the report, it would make sense to either just print a report with nicer formatting indicating where the errors are, or to return a data structure that lists the specific errors (not as a string, e.g. a list of dict(block=1, event='adc', field='delay', value=..., rounded_value=...)), and provide a function to nicely print that list.

[btasdelen]

I did not realize seq.check_timing() checks for ADC timing, good to know!

I realized that it also assumes rf_raster_time and adc_raster_time are the same?

pypulseq/pypulseq/check_timing.py

Lines 58 to 59 in 25938e8

	if hasattr(e, "type") and (e.type == "adc" or e.type == "rf"):
	raster = system.rf_raster_time

A dict is also awesome if one wants to easily search for certain events in a spam of errors.

[btasdelen]

I believe the answers and #200 addresses the issue. Closing.