[mcencini]

This PR adds standalone pypulseq.calc_rf_power() and the corresponding pypulseq.Sequence.calc_rf_power() method.

Compared to upstream MATLAB repository, this function by default directly uses explicitly the rf time vector to handle cases where pulse has not an uniform raster equal to rf raster (e.g., for block pulses). In addition, the Sequence method uses time_range instead of block_range to align with the other PyPulseq functions (e.g., plot().

This does not rely on any v1.5.1 features and can fit v1.5.0.

[github-actions]

Coverage Report

File	Stmts	Miss	Cover	Missing
/home/runner/.local/lib/python3.12/site-packages/pypulseq
add_gradients.py	137	21	85%	44, 52, 61, 85, 92, 125–128, 155, 162, 177, 187–190, 211–218, 225, 230, 261
add_ramps.py	36	36	0%	1–92
align.py	35	4	89%	41, 45, 69, 73
calc_duration.py	33	2	94%	43, 50
calc_ramp.py	220	216	2%	48–359
calc_rf_bandwidth.py	37	28	24%	45–81, 85–89
check_timing.py	99	28	72%	78, 82, 184, 203, 236, 243, 253–297
compress_shape.py	30	1	97%	28
convert.py	40	8	80%	42, 48, 66, 72–73, 82, 88–89
event_lib.py	96	15	84%	6–9, 48–51, 70–71, 185, 205–210
make_adc.py	108	10	91%	114, 163, 166, 170, 176, 180, 219, 221, 223, 231
make_adiabatic_pulse.py	150	39	74%	214–218, 238–242, 250–251, 274, 280, 349–368, 472–481, 519–527
make_arbitrary_grad.py	62	2	97%	124, 145
make_arbitrary_rf.py	83	24	71%	145–149, 152–156, 166–185, 188, 191
make_block_pulse.py	53	1	98%	132
make_delay.py	9	1	89%	27
make_extended_trapezoid.py	64	6	91%	74, 77, 137, 154, 156, 159
make_extended_trapezoid_area.py	110	4	96%	71, 74, 273, 276
make_gauss_pulse.py	90	5	94%	193, 201–203, 218
make_sigpy_pulse.py	129	28	78%	12–13, 134, 137, 186–190, 197–198, 201–202, 217, 224, 229, 267, 270, 295–305, 345, 348, 378–388
make_sinc_pulse.py	84	4	95%	196–198, 213
make_soft_delay.py	32	1	97%	136
make_trapezoid.py	128	7	95%	207, 221, 227, 245, 263, 268, 286
opts.py	66	9	86%	78, 83, 102, 142, 166–170
points_to_waveform.py	9	1	89%	27
rotate.py	69	5	93%	55, 69, 112, 119–120
rotate_3d.py	78	2	97%	145–146
sigpy_pulse_opts.py	28	7	75%	66–73
split_gradient.py	39	31	21%	46–103
split_gradient_at.py	70	27	61%	63–90, 110, 114, 118–120, 154–156
traj_to_grad.py	13	9	31%	26–40
/home/runner/.local/lib/python3.12/site-packages/pypulseq/SAR
SAR_calc.py	5	1	80%	9
/home/runner/.local/lib/python3.12/site-packages/pypulseq/Sequence
block.py	452	54	88%	64, 67, 75, 81, 96, 104, 110, 121, 124, 127, 135, 140, 149, 160, 168, 208, 210, 214, 329, 333, 336, 357, 372, 382, 393, 402, 434, 438, 447, 460, 471, 480–484, 502, 534–541, 558, 568, 594, 632, 671, 685, 712, 791, 828, 852
calc_grad_spectrum.py	81	76	6%	68–190
calc_pns.py	40	31	22%	45–96
ext_grad_check.py	55	2	96%	51, 53
ext_test_report.py	167	11	93%	79, 156, 167–168, 337–343
install.py	75	42	44%	31, 52, 69, 71, 112–131, 148, 181–184, 200–212, 254–278
parula.py	4	2	50%	19–86
read_seq.py	408	41	90%	44–45, 110, 117, 130, 133–134, 138, 184, 223–224, 242, 415, 436–453, 516, 519, 604, 628, 668, 699, 715–719, 726, 837, 848
sequence.py	711	133	81%	10–13, 108–118, 139–152, 206, 271–274, 321, 365, 441, 468–473, 671, 699–702, 793, 821–830, 842, 864, 905–908, 937, 963–969, 1004, 1015–1016, 1022, 1033, 1039, 1041, 1049, 1082–1090, 1222, 1320, 1326, 1329, 1332, 1371, 1496–1509, 1567, 1589–1591, 1612, 1675, 1683, 1781, 1792–1805, 1874–1875, 1886–1904, 1928, 1958–1966, 2002, 2016–2026
write_seq.py	362	7	98%	45, 77–79, 313, 347–349, 499
/home/runner/.local/lib/python3.12/site-packages/pypulseq/utils
cumsum.py	14	1	93%	17
paper_plot.py	63	58	8%	49–132
safe_pns_prediction.py	126	113	10%	50–87, 102–189, 197–214, 222, 244–250, 279–286, 310–336, 344–383, 396–411, 415
seq_plot.py	407	186	54%	22, 106, 138–139, 159, 176–178, 183–225, 235–242, 249–314 Warning: Comment truncated due to GitHub's 65536 character limit

[schuenke]

hey @mcencini

I started reviewing this PR as well and noticed a few bugs already that I will fix.

Furthermore, there are a few things to discuss:

1. time_range vs. block_range: I like the idea, but IMO we should use the same API as in Matlab by default. Therefore, my suggestion would be to have a similar implementation as Matlab / PyPulseq-matlab-like, but in addition allow to specify time_range instead of block_range. Would this be fine for you?
2. dt=None vs. dt=1e-6 vs. dt=system.rf_raster_time: Same as above, I like the idea of allowing dt=None, but would suggest to use the same default settings as in Matlab. However, hard-coding dt=1e-6 somehow feels wrong to me as well and I would suggest to adjust this in Matlab and Python to use system.rf_raster_time instead. Does this make sense to you?

I will probably start implementing it locally as suggested above, but it would be great to get your feedback on this as soon as possible.

[mcencini]

hi @schuenke, thanks for looking into this!

1. I agree in general, but I suppose we need to collect similar discrepancies in the rest of codebase like @ckolbPTB also suggested. For instance, I think Sequence.plot() also accepts time_range only, and I think the historical reason was that in principle there is no constraint for blocks to be added in order, so we should be a bit careful.

2. yes, it totally makes sense. I think I also initially thought to replace dt with system directly, but decided to keep dt for consistency with upstream MATLAB codebase. Perhaps it would make sense to also "fix" this in the MATLAB toolbox (@m-a-x-i-m-z)?

[schuenke]

.
1. I agree in general, but I suppose we need to collect similar discrepancies in the rest of codebase like @ckolbPTB also suggested. For instance, I think Sequence.plot() also accepts time_range only, and I think the historical reason was that in principle there is no constraint for blocks to be added in order, so we should be a bit careful.

Yeah, we should definitely do that for the existing code base, but for new implementations (like this one), I would prefer to adjust to Matlab already and, if possible, allow PyPulseq-like arguments (e.g. time_range) optionally. Long term the goal is to bring Matlab and Python closer together and if we introduce new functionality that deviates from Matlab, we have to adjust it at some point in the future leading to the same issues we currently have with seq.plot(), rotate(), make_trapezoid() etc...

[mcencini]

agreed - let's use block_range then :D

[schuenke]

I replaced the np.trapezoid() call in 1b0c681 because this requires numpy>=2.0, which requires python>=3.10 ?!, which would be a major change... its still a one liner

[schuenke]

This was also tested against MATLAB for various cases (single block/sinc/gauss pulses. different dt, full sequences, different block_ranges, different window_durations. Everything matching perfectly.

A comparison against MATLAB using different time_range values is not possible, because this is currently a Python-only feature.

[mcencini]

I replaced the np.trapezoid() call in 1b0c681 because this requires numpy>=2.0, which requires python>=3.10 ?!, which would be a major change... its still a one liner

That's ok - IMO at some point we should let end-of-life python versions and numpy 1.x go, but that probably deserve some long-term planning

[schuenke]

@mcencini: Same as for #341 , it would be great if you would review my changes again.

@sairamgeethanath @btasdelen @ckolbPTB @fzimmermann89 maybe one of you also wanna have a quick look at this? I tried to provide helpful infos in the comments.

[mcencini]

@schuenke it is ok to be merged to me, unless someone else wants to review it again as well.

[m-a-x-i-m-z]

Hi!

concerning one of the comments above:

yes, it totally makes sense. I think I also initially thought to replace dt with system directly, but decided to keep dt for consistency with upstream MATLAB codebase. Perhaps it would make sense to also "fix" this in the MATLAB toolbox (@m-a-x-i-m-z)?

The challenge in the Matlab's code that the function mr.calcRfPower() belongs to the static object mr, which has no data properties and therefore has no access to system.rfRasterTime.

The only reasonable fix in this direction is to make sure that the function Sequence.calcRfPower() explicitly passes its system.rfRasterTime to each call off mr.calcRfPower().

Should we do it this way?

--M.