A few small bugfixes

pypulseq/Sequence/block.py

@@ -146,6 +146,9 @@ def set_block(self, block_index: int, *args: SimpleNamespace) -> None:
                     "ref": label_id,
                 }
                 extensions.append(ext)
+        else:
+            # Floating point number given as delay
+            duration = max(duration, event)
 
     # =========
     # ADD EXTENSIONS

pypulseq/Sequence/calc_pns.py

@@ -53,7 +53,7 @@ def calc_pns(
     # 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)
+    max_t = max(g.x[-1] for g in gw_pp if g != None) - 1e-10
 
     # Determine sampling points
     if time_range == None:

pypulseq/Sequence/read_seq.py

@@ -384,7 +384,7 @@ def __read_definitions(input_file) -> Dict[str, str]:
                 value = value[0]
             definitions[tok[0]] = value
         except ValueError:  # Try clause did not work!
-            definitions[tok[0]] = tok[1:]
+            definitions[tok[0]] = line[len(tok[0])+1:].strip()
         line = __strip_line(input_file)
 
     return definitions
@@ -535,7 +535,7 @@ def __read_and_parse_events(input_file, *args: callable) -> EventLibrary:
     line = __strip_line(input_file)
 
     while line != "" and line != "#":
-        datas = re.split("(\s+)", line)
+        datas = re.split(r"(\s+)", line)
         datas = [d for d in datas if d != " "]
         data = np.zeros(len(datas) - 1, dtype=np.int32)
         event_id = int(datas[0])

pypulseq/Sequence/sequence.py

@@ -265,8 +265,12 @@ def calculate_kspace(
             # "Sample" ramps for display purposes otherwise piecewise-linear display (plot) fails
             ii = np.flatnonzero(np.abs(gm_pp[i].c[0, :]) > eps)
 
-            starts = np.int64(np.floor(gm_pp[i].x[ii] / self.grad_raster_time))
-            ends = np.int64(np.ceil(gm_pp[i].x[ii+1] / self.grad_raster_time))
+            # Do nothing if there are no ramps
+            if ii.shape[0] == 0:
+                continue
+
+            starts = np.int64(np.floor((gm_pp[i].x[ii] + eps) / self.grad_raster_time))
+            ends = np.int64(np.ceil((gm_pp[i].x[ii+1] - eps) / self.grad_raster_time))
 
             # Create all ranges starts[0]:ends[0], starts[1]:ends[1], etc.
             lengths = ends-starts+1
@@ -278,7 +282,8 @@ def calculate_kspace(
             inds[start_inds] = np.concatenate(([starts[0]], np.diff(starts) - lengths[:-1] + 1))
 
             tc.append(np.cumsum(inds) * self.grad_raster_time)
-        tc = np.concatenate(tc)
+        if tc != []:
+            tc = np.concatenate(tc)
 
 
         t_acc = 1e-10  # Temporal accuracy
@@ -588,7 +593,7 @@ def evaluate_labels(
         # Convert evolutions into label dictionary
         if len(label_evolution) > 0:
             for lab in labels:
-                labels[lab] = np.array([e[lab] for e in label_evolution])
+                labels[lab] = np.array([e[lab] if lab in e else 0 for e in label_evolution])
 
         return labels
 
@@ -773,15 +778,15 @@ def get_gradients(self,
                 raise Warning("Not all elements of the generated waveform are finite.")
 
             teps = 1e-12
-            if gw[0, 0] > 0 and gw[0, -1] < total_duration:
+            if gw[0, 0] > 0 and gw[0, -1] < total_duration - teps:
                 # teps terms to avoid integration errors over extended periods of time
                 _temp1 = np.array(([-teps, gw[0, 0] - teps], [0, 0]))
                 _temp2 = np.array(([gw[0, -1] + teps, total_duration + teps], [0, 0]))
                 gw = np.hstack((_temp1, gw, _temp2))
-            elif gw[0, 0] > 0:
+            elif gw[0, 0] > teps:
                 _temp = np.array(([-teps, gw[0, 0] - teps], [0, 0]))
                 gw = np.hstack((_temp, gw))
-            elif gw[0, -1] < total_duration:
+            elif gw[0, -1] < total_duration - teps:
                 _temp = np.array(([gw[0, -1] + teps, total_duration + teps], [0, 0]))
                 gw = np.hstack((gw, _temp))
 
@@ -1086,7 +1091,7 @@ def read(self, file_path: str, detect_rf_use: bool = False) -> None:
         read(self, path=file_path, detect_rf_use=detect_rf_use)
 
         # Initialize next free block ID
-        self.next_free_block_ID = max(self.block_events) + 1
+        self.next_free_block_ID = (max(self.block_events) + 1) if self.block_events else 1
 
     def register_adc_event(self, event: EventLibrary) -> int:
         return block.register_adc_event(self, event)

pypulseq/make_sigpy_pulse.py

@@ -114,6 +114,7 @@ def sigpy_n_seq(
     rfp.type = "rf"
     rfp.signal = signal
     rfp.t = t
+    rfp.shape_dur = t[-1]
     rfp.freq_offset = freq_offset
     rfp.phase_offset = phase_offset
     rfp.dead_time = system.rf_dead_time

pypulseq/make_trapezoid.py

@@ -42,7 +42,7 @@ def make_trapezoid(
     delay: float = 0,
     duration: float = 0,
     fall_time: float = 0,
-    flat_area: float = 0,
+    flat_area: float = None,
     flat_time: float = -1,
     max_grad: float = 0,
     max_slew: float = 0,
@@ -129,7 +129,7 @@ def make_trapezoid(
     else:
         fall_time = 0.0
 
-    if area is None and flat_area == 0 and amplitude == 0:
+    if area is None and flat_area is None and amplitude == 0:
         raise ValueError("Must supply either 'area', 'flat_area' or 'amplitude'.")
 
     if flat_time != -1:
@@ -139,7 +139,7 @@ def make_trapezoid(
                 and rise_time > 0:
             # We have rise_time, flat_time and area.
             amplitude2 = area / (rise_time + flat_time)
-        elif flat_area > 0:
+        elif flat_area is not None:
             amplitude2 = flat_area / flat_time
         else:
             raise ValueError(

pypulseq/tests/test_sigpy.py

@@ -6,6 +6,7 @@
 import numpy as np
 import sigpy.mri.rf as rf
 
+import pypulseq as pp
 from pypulseq.make_sigpy_pulse import sigpy_n_seq
 from pypulseq.opts import Opts
 from pypulseq.sigpy_pulse_opts import SigpyPulseOpts
@@ -48,6 +49,9 @@ def test_slr(self):
             pulse_cfg=pulse_cfg,
             plot=False,
         )
+
+        seq = pp.Sequence()
+        seq.add_block(rfp)
 
         [a, b] = rf.sim.abrm(
             pulse,
@@ -100,6 +104,9 @@ def test_sms(self):
             pulse_cfg=pulse_cfg,
             plot=False
         )
+
+        seq = pp.Sequence()
+        seq.add_block(rfp)
 
         [a, b] = rf.sim.abrm(
             pulse,

pypulseq/traj_to_grad.py

@@ -29,14 +29,14 @@ def traj_to_grad(
         raster_time = Opts.default.grad_raster_time
 
     # Compute finite difference for gradients in Hz/m
-    g = (k[1:] - k[:-1]) / raster_time
+    g = (k[...,1:] - k[...,:-1]) / raster_time
     # Compute the slew rate
-    sr0 = (g[1:] - g[:-1]) / raster_time
+    sr0 = (g[...,1:] - g[...,:-1]) / raster_time
 
     # Gradient is now sampled between k-space points whilst the slew rate is between gradient points
-    sr = np.zeros(len(sr0) + 1)
-    sr[0] = sr0[0]
-    sr[1:-1] = 0.5 * (sr0[-1] + sr0[1:])
-    sr[-1] = sr0[-1]
+    sr = np.zeros(sr0.shape[:-1] + (sr0.shape[-1] + 1,))
+    sr[...,0] = sr0[...,0]
+    sr[...,1:-1] = 0.5 * (sr0[...,:-1] + sr0[...,1:])
+    sr[...,-1] = sr0[...,-1]
 
     return g, sr

pypulseq/utils/siemens/readasc.py

@@ -40,7 +40,7 @@ def readasc(filename : str) -> Tuple[dict, dict]:
             # regex wizardry: Matches lines like 'a[0].b[2][3].c = "string" # comment'
             # Note this assumes correct formatting, e.g. does not check whether
             # brackets match.
-            match = re.match('^\s*([a-zA-Z0-9\[\]\._]+)\s*\=\s*((\"[^\"]*\"|\\\'[^\\\']\\\')|(\d+)|([0-9\.e\-]+))\s*((#|\/\/)(.*))?$', next_line)
+            match = re.match(r'^\s*([a-zA-Z0-9\[\]\._]+)\s*\=\s*(("[^"]*"|\'[^\']\')|(\d+)|([0-9\.e\-]+))\s*((#|\/\/)(.*))?$', next_line)
 
             if match:
                 field_name = match[1]