Create your own sequence PyPulseq version - josalggui/MaRGE GitHub Wiki

Adding a New Sequence to the GUI in MaRGE

Starting with version v0.6.0, MaRGE allows users to implement custom MRI sequences using PyPulseq. Writing sequences compatible with MaRGE follows a structure similar to the mriBlankSeq framework. All sequences still inherit from mriBlankSeq, which now provides methods to convert PyPulseq sequences into mriBlankSeq arrays.

Sequence script must be saved in marge/seq folder. A ready-to-use sequence template is available at marge/seq/seq_template.py

This template includes instructions for making your sequence compatible with MaRGE, allowing it to be executed directly from the graphical user interface (GUI).

The overall structure remains the same as in the standard mriBlankSeq framework. Each sequence must implement the following methods:

  • sequenceInfo
  • sequenceTime
  • sequenceAtributes
  • sequenceRun
  • sequenceAnalysis (optional if a standalone analysis is created)

The primary difference is in the sequenceRun method, where the sequence is implemented using PyPulseq.

Step 1: Define the interpreter for FloSeq/PSInterpreter

The interpreter converts the high-level pulse sequence description into low-level instructions that can be executed by the scanner hardware. This interpreter is internally configured during scanner calibration.

self.flo_interpreter = PSInterpreter(
    tx_warmup=hw.blkTime,  # Transmit chain warm-up time (us)
    rf_center=hw.larmorFreq * 1e6,  # Larmor frequency (Hz)
    rf_amp_max=hw.b1Efficiency / (2 * np.pi) * 1e6,  # Maximum RF amplitude (Hz)
    gx_max=hw.gFactor[0] * hw.gammaB,  # Maximum gradient amplitude in X (Hz/m)
    gy_max=hw.gFactor[1] * hw.gammaB,  # Maximum gradient amplitude in Y (Hz/m)
    gz_max=hw.gFactor[2] * hw.gammaB,  # Maximum gradient amplitude in Z (Hz/m)
    grad_max=np.max(hw.gFactor) * hw.gammaB,  # Maximum gradient amplitude (Hz/m)
    grad_t=hw.grad_raster_time * 1e6  # Gradient raster time (us)
)

Step 2: Define system properties using PyPulseq (pp.Opts)

These parameters describe the hardware capabilities of the MRI system, including gradient limits, slew rates, and timing constraints. They are typically derived from the hardware configuration.

self.system = pp.Opts(
    rf_dead_time=hw.blkTime * 1e-6,  # Dead time between RF pulses (s)
    max_grad=np.max(hw.gFactor) * 1e3,  # Maximum gradient strength (mT/m)
    grad_unit='mT/m',  # Units of gradient strength
    max_slew=hw.max_slew_rate,  # Maximum gradient slew rate (mT/m/ms)
    slew_unit='mT/m/ms',  # Units of gradient slew rate
    grad_raster_time=hw.grad_raster_time,  # Gradient raster time (s)
    rise_time=hw.grad_rise_time,  # Gradient rise time (s)
    rf_raster_time=1e-6,
    block_duration_raster=1e-6
)

Step 3: Perform any necessary calculations for the sequence

Before defining the sequence blocks, compute all necessary parameters, such as:

  • Timing delays
  • RF amplitudes
  • Gradient strengths
  • Readout durations

These values will be used when constructing the PyPulseq blocks.

Step 4: Define the experiment to obtain the true bandwidth

To determine the actual acquisition bandwidth, initialize an experiment and retrieve the sampling period using get_sampling_period().

import marge.controller.experiment_gui as ex
if not self.demo:
    expt = ex.Experiment(
        lo_freq=hw.larmorFreq,  # Larmor frequency in MHz
        rx_t=sampling_period,  # Sampling time in us
        init_gpa=False,  # Whether to initialize GPA board (False for True)
        gpa_fhdo_offset_time=(1 / 0.2 / 3.1),  # GPA offset time calculation
        auto_leds=True  # Automatic control of LEDs (False or True)
    )
    sampling_period = expt.get_sampling_period()  # us
    bw = 1 / sampling_period  # MHz
    print("Acquisition bandwidth fixed to: %0.3f kHz" % (bw * 1e3))
    expt.__del__()
self.mapVals['bw_MHz'] = bw
self.mapVals['sampling_period_us'] = sampling_period

More information about experiment_gui class here.

Step 5: Define the sequence blocks

Define the fundamental building blocks of the MRI sequence, including:

  • RF pulses
  • Gradient waveforms
  • ADC events

Note
The conversion from PyPulseq to mriBlankSeq automatically applies the gradient delay specified in hardware specifications. You do not need to manually compensate for gradient delays. However, ensure that the first gradient event does not begin earlier than the configured gradient delay, as this would result in negative timing and cause a sequence error.

Step 6: Define initializeBatch

In this step, implement the initializeBatch method, which is responsible for creating and initializing a new batch of sequence blocks. This method is typically used to insert dummy pulses that are executed at the beginning of each batch to ensure steady-state conditions.

The method should return:

  • A PyPulseq Sequence object containing the initialized blocks
  • A dictionary containing the total number of acquired points for each batch.
  • The total number of ADC events in the sequence

Example Implementation

def initializeBatch():
    """
    Initializes a batch of MRI sequence blocks using PyPulseq for a given
    experimental configuration.

    Returns
    -------
    tuple
        - batch (pp.Sequence): PyPulseq sequence object containing the batch.
        - n_rd_points (int): Total number of readout points in the batch.
        - n_adc (int): Total number of ADC acquisitions in the batch.
    """

    # Instantiate PyPulseq sequence object
    batch = pp.Sequence(system)
    n_rd_points = 0
    n_adc = 0

    # Add dummy pulses to the batch
    for _ in range(self.dummyPulses):
        batch.add_block(rf_ex, delay_repetition)

    return batch, n_rd_points, n_adc

Each time a new batch is created, this method is called to ensure that the required dummy pulses are included before any data-acquisition blocks are added.

Step 7: Define createBatches

In this step, implement the createBatches method. This method is responsible for assembling the full pulse sequence by populating batches with the sequence blocks defined earlier and managing hardware constraints on the maximum number of acquired readout points.

The sequence is split into multiple batches when the total number of readout points exceeds the hardware limit (hw.maxRdPoints). Each batch is independently written, interpreted, and executed.

Responsibilities of createBatches

  • Iterate over all repetitions of the sequence (e.g., slices or phase encodes)
  • Add RF, ADC, and delay blocks to the active batch
  • Track the total number of acquired readout points
  • Start a new batch when the hardware readout limit is exceeded
  • Write and interpret each batch using the FloSeq interpreter

Example Implementation

def createBatches():
    """
    Create batches for the full pulse sequence.

    Instructions
    ------------
    - This function builds the complete pulse sequence by iterating over
      all repetitions.
    - Each iteration adds the RF pulse, ADC block, and repetition delay.
    - When the maximum number of readout points is exceeded, a new batch
      is created and the previous one is finalized.

    Returns
    -------
    waveforms : dict
        Dictionary containing the generated waveforms for each batch.
    n_rd_points_dict : dict
        Number of readout points per batch.
    n_adc : int
        Total number of ADC acquisitions across all batches.
    """

    batches = {}            # PyPulseq sequence objects per batch
    waveforms = {}          # Generated waveforms per batch
    n_rd_points_dict = {}   # Readout points per batch

    n_rd_points = 0         # Readout points counter for current batch
    seq_idx = 0             # Batch index
    n_adc = 0               # Total ADC counter
    batch_num = "batch_0"   # Initial batch name

    # Loop over all repetitions
    for repetition in range(self.nRepetitions):

        # Check if a new batch is required
        if seq_idx == 0 or n_rd_points + self.nPoints > hw.maxRdPoints:

            # Finalize previous batch if it exists
            if seq_idx > 0:
                batches[batch_num].write(batch_num + ".seq")
                waveforms[batch_num], param_dict = flo_interpreter.interpret(
                    batch_num + ".seq"
                )
                print(f"{batch_num}.seq ready!")

            # Initialize a new batch
            seq_idx += 1
            n_rd_points_dict[batch_num] = n_rd_points
            n_rd_points = 0
            batch_num = f"batch_{seq_idx}"

            batches[batch_num], n_rd_points, n_adc_0 = initializeBatch()
            n_adc += n_adc_0
            print(f"Creating {batch_num}.seq...")

        # Add sequence blocks to the current batch
        batches[batch_num].add_block(rf_ex, adc, delay_repetition)
        n_rd_points += self.nPoints
        n_adc += 1

    # Finalize the last batch
    batches[batch_num].write(batch_num + ".seq")
    waveforms[batch_num], param_dict = flo_interpreter.interpret(batch_num + ".seq")
    print(f"{batch_num}.seq ready!")
    print(f"{len(batches)} batches created. Sequence ready!")

    # Update readout points for the final batch
    n_rd_points_dict.pop('batch_0')
    n_rd_points_dict[batch_num] = n_rd_points

    return waveforms, n_rd_points_dict, n_adc

This approach ensures that long sequences are safely split into hardware-compatible batches while preserving correct timing, acquisition counts, and waveform generation.

Step 8: Run the batches

This step executes the batches and retrieves the acquired data. The oversampled data will be available in self.mapVals['data_over'], and the decimated data in self.mapVals['data_decimated

waveforms, n_readouts, n_adc = createBatches()
return self.runBatches(waveforms=waveforms,
                       n_readouts=n_readouts,
                       n_adc=n_adc,
                       frequency=hw.larmorFreq,  # MHz
                       bandwidth=bw,  # MHz
                       decimate='Normal',
                       )

More info about runBatches method can be found here.