Merge pull request #159 from FrankZijlstra/bugfixes

Various bug fixes, features, and cleanups

pypulseq/Sequence/block.py

@@ -97,17 +97,6 @@ def set_block(self, block_index: int, *args: SimpleNamespace) -> None:
                 if self.block_events[block_index][idx] != 0:
                     raise ValueError(f'Multiple {event.channel.upper()} gradient events were specified in set_block')
 
-                check_g[channel_num] = SimpleNamespace()
-                check_g[channel_num].idx = idx
-                check_g[channel_num].start = (0, 0)
-                check_g[channel_num].stop = (
-                        event.delay
-                        + event.rise_time
-                        + event.fall_time
-                        + event.flat_time,
-                        0,
-                    )
-
                 if hasattr(event, "id"):
                     trap_id = event.id
                 else:
@@ -206,6 +195,9 @@ def set_block(self, block_index: int, *args: SimpleNamespace) -> None:
                 )
 
             if block_index > 1:
+                # TODO: block_index - 1 assumes linear ordering of blocks.
+                #       However, searching for the previous block index in block_events is a big performance hit
+                #       For newly inserted blocks (block_index not in self.block_events), we can just use the last block as previous
                 prev_id = self.block_events[block_index - 1][grad_to_check.idx]
                 if prev_id != 0:
                     prev_lib = self.grad_library.get(prev_id)
@@ -281,13 +273,13 @@ def get_block(self, block_index: int) -> SimpleNamespace:
         block.delay = delay
 
     if event_ind[1] > 0:  # RF
-        if len(self.rf_library.type) >= event_ind[1]:
+        if event_ind[1] in self.rf_library.type:
             block.rf = self.rf_from_lib_data(
                 self.rf_library.data[event_ind[1]], self.rf_library.type[event_ind[1]]
             )
         else:
             block.rf = self.rf_from_lib_data(
-                self.rf_library.data[event_ind[1]]
+                self.rf_library.data[event_ind[1]] , 'u'
             )  # Undefined type/use
 
     # Gradients

pypulseq/Sequence/calc_pns.py

@@ -1,6 +1,6 @@
 import math
 from types import SimpleNamespace
-from typing import Tuple
+from typing import Tuple, List
 
 import matplotlib.pyplot as plt
 import pypulseq as pp
@@ -14,7 +14,10 @@
 
 
 def calc_pns(
-        obj : Sequence, hardware : SimpleNamespace, do_plots: bool = True
+        obj: Sequence,
+        hardware: SimpleNamespace,
+        time_range: List[float] = None,
+        do_plots: bool = True
         ) -> Tuple[bool, np.array, np.ndarray, np.array]:
     """
     Calculate PNS using safe model implementation by Szczepankiewicz and Witzel
@@ -46,44 +49,42 @@ def calc_pns(
         Time axis for the pns_norm and pns_components arrays
     """
 
-    # acquire the entire gradient wave form
-    gw = obj.waveforms_and_times()[0]
+    dt = obj.grad_raster_time
+    # Get gradients as piecewise-polynomials
+    gw_pp = obj.get_gradients(time_range=time_range)
+    ng = len(gw_pp)
+    max_t = max(g.x[-1] for g in gw_pp if g != None)
+
+    # Determine sampling points
+    if time_range == None:
+        nt = int(np.ceil(max_t/dt))
+        t = (np.arange(nt) + 0.5)*dt
+    else:
+        tmax = min(time_range[1], max_t) - max(time_range[0], 0)
+        nt = int(np.ceil(tmax/dt))
+        t = max(time_range[0], 0) + (np.arange(nt) + 0.5)*dt
+
+    # Sample gradients
+    gw = np.zeros((t.shape[0], ng))
+    for i in range(ng):
+        if gw_pp[i] != None:
+            gw[:,i] = gw_pp[i](t)
+
+
     if do_plots:
         plt.figure()
-        plt.plot(gw[0][0], gw[0][1], gw[1][0], gw[1][1], gw[2][0], gw[2][1]) # plot the entire gradient shape
-        plt.title('gradient wave form, in T/m')
-
-    # find beginning and end times and resample GWs to a regular sampling raster
-    tf = []
-    tl = []
-    for i in range(3):
-        if gw[i].shape[1] > 0:
-            tf.append(gw[i][0,0])
-            tl.append(gw[i][0,-1])
-
-    nt_min = math.floor(min(tf) / obj.grad_raster_time + pp.eps) 
-    nt_max = math.ceil(max(tl) / obj.grad_raster_time - pp.eps)
+        for i in range(ng):
+            if gw_pp[i] != None:
+                plt.plot(gw_pp[i].x[1:-1], gw_pp[i].c[1,:-1])
+        plt.title('gradient wave form, in Hz/m')
 
-    # shift raster positions to the centers of the raster periods
-    nt_min = nt_min + 0.5
-    nt_max = nt_max - 0.5
-    if nt_min < 0.5:
-        nt_min = 0.5
-
-    t_axis = (np.arange(0, math.floor(nt_max-nt_min) + 1) + nt_min) * obj.grad_raster_time
-
-    gwr = np.zeros((t_axis.shape[0],3))
-    for i in range(3):
-        if gw[i].shape[1] > 0:
-            gwr[:,i] = np.interp(t_axis, gw[i][0], gw[i][1])
-
     if type(hardware) == str:
         # this loads the parameters from the provided text file
         asc, _ = readasc(hardware)
         hardware = asc_to_hw(asc)
 
     # use the Szczepankiewicz' and Witzel's implementation
-    [pns_comp,res] = safe_gwf_to_pns(gwr/obj.system.gamma, np.nan*np.ones(t_axis.shape[0]), obj.grad_raster_time, hardware) # the RF vector is unused in the code inside but it is zeropaded and exported ... 
+    [pns_comp,res] = safe_gwf_to_pns(gw/obj.system.gamma, np.nan*np.ones(t.shape[0]), obj.grad_raster_time, hardware) # the RF vector is unused in the code inside but it is zeropaded and exported ... 
 
     # use the exported RF vector to detect and undo zero-padding
     pns_comp = 0.01 * pns_comp[~np.isfinite(res.rf[1:]),:]
@@ -98,4 +99,4 @@ def calc_pns(
         plt.figure()
         safe_plot(pns_comp*100, obj.grad_raster_time)
 
-    return ok, pns_norm, pns_comp, t_axis
+    return ok, pns_norm, pns_comp, t

pypulseq/Sequence/ext_test_report.py

@@ -15,7 +15,7 @@ def ext_test_report(self) -> str:
     """
     # Find RF pulses and list flip angles
     flip_angles_deg = []
-    for k in self.rf_library.keys:
+    for k in self.rf_library.data:
         lib_data = self.rf_library.data[k]
         if len(self.rf_library.type) >= k:
             rf = self.rf_from_lib_data(lib_data, self.rf_library.type[k])

pypulseq/Sequence/read_seq.py

@@ -215,7 +215,7 @@ def read(self, path: str, detect_rf_use: bool = False) -> None:
     # Fix blocks, gradients and RF objects imported from older versions
     if version_combined < 1004000:
         # Scan through RF objects
-        for i in self.rf_library.data.keys():
+        for i in self.rf_library.data:
             self.rf_library.data[i] = [
                 *self.rf_library.data[i][:3],
                 0,
@@ -224,7 +224,7 @@ def read(self, path: str, detect_rf_use: bool = False) -> None:
             self.rf_library.lengths[i] += 1
 
         # Scan through the gradient objects and update 't'-s (trapezoids) und 'g'-s (free-shape gradients)
-        for i in self.grad_library.data.keys():
+        for i in self.grad_library.data:
             if self.grad_library.type[i] == "t":
                 if self.grad_library.data[i][1] == 0:
                     if (
@@ -253,24 +253,25 @@ def read(self, path: str, detect_rf_use: bool = False) -> None:
         # For versions prior to 1.4.0 block_durations have not been initialized
         self.block_durations = dict()
         # Scan through blocks and calculate durations
-        for block_counter in range(len(self.block_events)):
-            block = self.get_block(block_counter + 1)
-            if delay_ind_temp[block_counter + 1] > 0:
+        for block_counter in self.block_events:
+            block = self.get_block(block_counter)
+            if delay_ind_temp[block_counter] > 0:
                 block.delay = SimpleNamespace()
                 block.delay.type = "delay"
                 block.delay.delay = temp_delay_library.data[
-                    delay_ind_temp[block_counter + 1]
+                    delay_ind_temp[block_counter]
                 ]
-            self.block_durations[block_counter + 1] = calc_duration(block)
+            self.block_durations[block_counter] = calc_duration(block)
 
+    # TODO: Is it possible to avoid expensive get_block calls here? 
     grad_channels = ["gx", "gy", "gz"]
     grad_prev_last = np.zeros(len(grad_channels))
-    for block_counter in range(len(self.block_events)):
-        block = self.get_block(block_counter + 1)
+    for block_counter in self.block_events:
+        block = self.get_block(block_counter)
         block_duration = block.block_duration
         # We also need to keep track of the event IDs because some PyPulseq files written by external software may contain
         # repeated entries so searching by content will fail
-        event_idx = self.block_events[block_counter + 1]
+        event_idx = self.block_events[block_counter]
         # Update the objects by filling in the fields not contained in the PyPulseq file
         for j in range(len(grad_channels)):
             grad = getattr(block, grad_channels[j])
@@ -283,16 +284,15 @@ def read(self, path: str, detect_rf_use: bool = False) -> None:
                     grad_prev_last[j] = 0
 
                 if hasattr(grad, "first"):
+                    grad_prev_last[j] = grad.last
                     continue
 
                 amplitude_ID = event_idx[j + 2]
                 if amplitude_ID in event_idx[2:(j+2)]: # We did this update for the previous channels, don't do it again.
-
-                    # Remove the block from the cache and re-add it, otherwise, the other grad with same id in the block will not be updated.
                     if self.use_block_cache:
-                        del self.block_cache[block_counter + 1] # TODO: I feel like the update_data function should be responsible for this.
-                        block = self.get_block(block_counter + 1)
-
+                        # Update block cache in-place using the first/last values that should now be in the grad_library
+                        grad.first = self.grad_library.data[amplitude_ID][4]
+                        grad.last = self.grad_library.data[amplitude_ID][5]
                     continue
 
                 grad.first = grad_prev_last[j]
@@ -301,6 +301,7 @@ def read(self, path: str, detect_rf_use: bool = False) -> None:
                     grad_duration = grad.delay + grad.tt[-1]
                 else:
                     # Restore samples on the edges of the gradient raster intervals for that we need the first sample
+                    # TODO: This code does not always restore reasonable values for grad.last
                     odd_step1 = [grad.first, *2 * grad.waveform]
                     odd_step2 = odd_step1 * (np.mod(range(len(odd_step1)), 2) * 2 - 1)
                     waveform_odd_rest = np.cumsum(odd_step2) * (
@@ -318,31 +319,23 @@ def read(self, path: str, detect_rf_use: bool = False) -> None:
                     grad_prev_last[j] = 0
 
                 amplitude = self.grad_library.data[amplitude_ID][0]
-                if version_combined >= 1004000:
-                    old_data = np.array(
-                        [amplitude, grad.shape_id, grad.time_id, grad.delay]
-                    )
-                else:
-                    old_data = np.array([amplitude, grad.shape_id, grad.delay])
-                new_data = np.array(
-                    [
+                new_data = (
                         amplitude,
                         grad.shape_id,
                         grad.time_id,
                         grad.delay,
                         grad.first,
                         grad.last,
-                    ]
                 )
-                self.grad_library.update_data(amplitude_ID, old_data, new_data, "g")
+                self.grad_library.update_data(amplitude_ID, None, new_data, "g")
 
             else:
                 grad_prev_last[j] = 0
 
     if detect_rf_use:
         # Find the RF pulses, list flip angles, and work around the current (rev 1.2.0) Pulseq file format limitation
         # that the RF pulse use is not stored in the file
-        for k in self.rf_library.keys.keys():
+        for k in self.rf_library.data:
             lib_data = self.rf_library.data[k]
             rf = self.rf_from_lib_data(lib_data)
             flip_deg = np.abs(np.sum(rf.signal[:-1] * (rf.t[1:] - rf.t[:-1]))) * 360
@@ -359,7 +352,12 @@ def read(self, path: str, detect_rf_use: bool = False) -> None:
                 else:
                     self.rf_library.type[k] = "r"
             self.rf_library.data[k] = lib_data
-
+
+            # Clear block cache for all blocks that contain the modified RF event
+            for block_counter,events in self.block_events.items():
+                if events[1] == k:
+                    del self.block_cache[block_counter]
+
 
 def __read_definitions(input_file) -> Dict[str, str]:
     """