Usage
=====

Running a Pulse Sequence
------------------------

Pulse sequence programs are used programmatically via the :class:`~matipo.sequence.Sequence` class. Below is the code for a typical usage scenario. This code may be run directly in a jupyter notebook.

Load a pulse sequence by passing its file path::

    from matipo import SEQUENCE_DIR, GLOBALS_DIR
    from matipo.sequence import Sequence

    # load system FID pulse sequence
    fid_seq = Sequence(SEQUENCE_DIR+'FID.py')

Load parameter files by passing their file paths to the :meth:`~matipo.sequence.Sequence.loadpar` method::

    # load system frequency and pulse calibration files
    fid_seq.loadpar(GLOBALS_DIR+'frequency.yaml')
    fid_seq.loadpar(GLOBALS_DIR+'hardpulse_90.yaml')

Set parameters directly using the :meth:`~matipo.sequence.Sequence.setpar` method::

    fid_seq.setpar(
        n_scans=1,
        n_samples=10000,
        t_dw=10e-6,
        t_end=1
    )

The current set of parameters can be read using the :attr:`~matipo.sequence.Sequence.par` property::

    # print all parameters
    print(fid_seq.par)
    # use parameters to calculate total acquisition time: (dwell time)*(samples)
    print('Acquisition time:', fid_seq.par.t_dw*fid_seq.par.n_samples)

.. note::
    All parameters have a default value set in the pulse sequence program which will be used if not overriden by :meth:`~matipo.sequence.Sequence.loadpar` or :meth:`~matipo.sequence.Sequence.setpar`.

Run the pulse sequence program using the :meth:`~matipo.sequence.Sequence.run` method, which also returns the data. This is an ``async`` method, which means that ``await`` must be used to schedule it to run and obtain the result when finished::

    data = await fid_seq.run()

The data from a sequence object's last run can also be accessed later with the :attr:`~matipo.sequence.Sequence.data` property::

    data = fid_seq.data

Data is returned as a numpy array of complex values. This data may be plotted with matplotlib::

    import matplotlib.pyplot as plt
    plt.plot(data.real, label='real')
    plt.plot(data.imag, label='imag')
    plt.legend()
    plt.show()

Parameters may be saved to a file using the :meth:`~matipo.sequence.Sequence.savepar` method::

    fid_seq.savepar('FID_example.yaml')

This file can then be loaded with :meth:`~matipo.sequence.Sequence.loadpar` to restore the same parameters at a later date.

Data may be saved using numpy::

    import numpy as np
    np.save('FID_example.npy', data)

Writing Custom Pulse Sequence Programs
--------------------------------------

Template
^^^^^^^^

There is some boilerplate that all sequence programs must have. Here is a simple sequence which runs a single acquisition on a channel, and can be used as a template::

    from matipo.sequence import ParDef
    from collections import namedtuple
    import numpy as np

    # Define input parameters (must be named 'PARDEF')
    PARDEF = [
        # args: name, dtype, default value
        ParDef('f', float, 1e6, unit='Hz'),
        ParDef('t_dw', float, 1e-6, min=0.1e-6, max=160e-6, unit='s'),
        ParDef('n_samples', int, 1000, min=2, unit=''),
        ParDef('rx_chan', int, 0, unit='')
    ]

    # Create parameter set type (must be named 'ParameterSet')
    ParameterSet = namedtuple('ParameterSet', [pd.name for pd in PARDEF])

    # Main pulse sequence function, taking the sequencer object and a parameter set
    def main(seq, par: ParameterSet):
        t_acq = par.n_samples * par.t_dw
        t_flatfilter_groupdelay = 10*par.t_dw

        yield seq.rx[par.rx_chan].freq(par.f)
        yield seq.rx[par.rx_chan].dwelltime(par.t_dw)
        yield seq.rx[par.rx_chan].mode(flatfilter=True)

        # Wait filter group delay after changing acquisition settings to avoid edge effects
        yield seq.wait(t_flatfilter_groupdelay)
        
        # Start acquisition with ID=0 and the desired number of samples
        yield seq.rx[par.rx_chan].acquire(0, par.n_samples)

        # Allow acquisition to finish completely before ending the sequence
        yield seq.wait(t_acq + t_flatfilter_groupdelay)

Input Parameter Definition
^^^^^^^^^^^^^^^^^^^^^^^^^^

The input parameters are defined in the :code:`PARDEF` list::

    PARDEF = [
        # args: name, dtype, default value
        ParDef('f', float, 1e6, unit='Hz'),
        ParDef('t_dw', float, 1e-6, min=0.1e-6, max=160e-6, unit='s'),
        ParDef('n_samples', int, 1000, min=2, unit=''),
        ParDef('rx_chan', int, 0, unit='')
    ]

Where each parameter is defined with :code:`ParDef(<name>, <data type>, <default value>)`. ``min``, ``max``, and ``unit`` may also be defined, and will be read by experiments using :class:`~matipo.experiment.BaseExperiment` to set limits on the inputs. The `Pint library <https://pint.readthedocs.io/>`_ is used to convert units if Pint Quantities are passed into the :meth:`~matipo.sequence.Sequence.setpar` method.

The :obj:`matipo.sequence.floatarray` data type may be used to pass numpy arrays as parameters.

Main Function
^^^^^^^^^^^^^

The :code:`main()` function is the actual pulse sequence program::

    def main(seq, par: ParameterSet):
        t_acq = par.n_samples * par.t_dw
        t_flatfilter_groupdelay = 10*par.t_dw

        yield seq.rx[par.rx_chan].freq(par.f)
        yield seq.rx[par.rx_chan].dwelltime(par.t_dw)
        yield seq.rx[par.rx_chan].mode(flatfilter=True)

        # Wait filter group delay after changing acquisition settings to avoid edge effects
        yield seq.wait(t_flatfilter_groupdelay)
        
        # Start acquisition with ID=0 and the desired number of samples
        yield seq.rx[par.rx_chan].acquire(0, par.n_samples)

        # Allow acquisition to finish completely before ending the sequence
        yield seq.wait(t_acq + t_flatfilter_groupdelay)

This function is a python generator which returns pulse sequence commands using Python's ``yield`` syntax.

Pulse Sequence Commands
^^^^^^^^^^^^^^^^^^^^^^^

Pulse sequence commands are methods of the Sequencer object and its children:

``seq`` (:class:`~matipo.hardware.sequence_inst.SequencerTop`):
    * ``wait(<time>)``
    * ``join(<list of commands>)``
    * ``tx[<channel index>]`` (:class:`~matipo.hardware.sequence_inst.RFTX`)
        * ``enable()``
        * ``disable()``
        * ``freq(<freq (Hz)>)``
        * ``phase(<phase (degrees)>)``
        * ``amp(<amplitude (ratio)>)``
    * ``rx[<channel index>]`` (:class:`~matipo.hardware.sequence_inst.ACQ`)
        * ``acquire(<id>, <number of samples>)``
        * ``dwelltime(<sample dwell time (s)>)``
        * ``freq(<freq (Hz)>)``
        * ``phase(<phase (degrees)>)``
        * ``mode(flatfilter={True|False}, raw={True|False})``
    * ``grad[<controller index>]`` (:class:`~matipo.hardware.sequence_inst.GRAD`)
        * ``vec(<x>, <y>, <z>)``
        * ``aux(<value>)``
    * ``shim[<controller index>]`` (:class:`~matipo.hardware.sequence_inst.SHIM`)
        * ``set(<channel>, <value>)``
    * ``gpo[<controller index>]`` (:class:`~matipo.hardware.sequence_inst.GPO`)
        * ``write(<value>)``
        * ``set(<bitmask>)``
        * ``clear(<bitmask>)``

Composing Commands using Variables & Functions
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A list of commands may be joined together with ``seq.join()`` and stored in a variable to be yielded multiple times for efficiency, e.g.::

    pulse_90 = seq.join([
        seq.tx[0].freq(10e6),
        seq.tx[0].amp(1),
        seq.tx[0].enable(),
        seq.wait(10e-6),
        seq.tx[0].disable(),
        seq.tx[0].amp(0)
    ])
    yield pulse_90
    yield seq.wait(10e-6)
    yield pulse_90
    yield seq.wait(10e-6)
    yield pulse_90
    yield seq.wait(10e-6)

This code will result in three 10 microsecond full amplitude pulses at 10 MHz with 10 microseconds spacing.

Custom functions may also be defined using python syntax for convenience, e.g. to transmit an RF pulse with amplitude modulation::

    def amplitude_modulated_pulse(seq, chan, amp_pts, duration):
        # calculate spacing of amplitude updates
        dt = duration/len(amp_pts)

        # start a list of commands with the TX enable command
        cmds = [seq.tx[chan].enable()]

        # for each amplitude value, append a command setting that amplitude and waiting for the spacing
        cmds += [
            seq.join([
                seq.tx[chan].amp(amp),
                seq.wait(dt)])
            for amp in amp_pts]

        # append TX disable command
        cmds += [seq.tx[chan].disable()]

        # join list of commands and return
        return seq.join(cmds)

    # example usage, 1 millisecond 5 lobe sinc pulse:
    def main(seq, par):
        ... snipped
        sinc_amp_pts = np.sinc(np.linspace(-3, 3, 100))
        yield amplitude_modulated_pulse(seq, 0, sinc_amp_pts, 1e-3)
        ... rest of sequence

Amplitude Ranges
^^^^^^^^^^^^^^^^
Transmit amplitude is a value in the range -1.0 to 1.0, where 1 is the maximum possible amplitude of 0dBm.

Gradient and shim amplitudes are also defined in the range -1.0 to 1.0, which corresponds to an output range of -10V to 10V for the gradient channels and -5V to 5V for the shim channels.

RF pulse amplitude may be negative for convenience, which results in a 180 degree phase shift. This makes defining shaped pulses simpler as the phase can be kept constant and the amplitude varied with a waveform that contains negative values.

Timing Constraints
^^^^^^^^^^^^^^^^^^

The minimum wait time is 100 nanoseconds (``100e-9``), but the timing precision is 10 nanoseconds. This means you may have a delay of exactly 110 nanoseconds, for example, but a wait time of 116 nanoseconds would be rounded to 120 nanoseconds.

There must be at least a :code:`10e-6` (10 microseconds) delay between successive updates on a single gradient/shim controller. Note that you may still have multiple commands before a time delay as they will be combined together into a single update to the controller by the driver.

If a timing contraint violation occurs, they driver will terminate the sequence and return an :ref:`error code <error-codes>`.

Acquisition Modes/Filters
^^^^^^^^^^^^^^^^^^^

The signal processing pipeline is a Digital-Down-Converter (DDC) with three stages:
    #. Frequency Mixer: Shifts the reference frame to the desired frequency and outputs quadrature data. The frequency and phase are set with the :meth:`~matipo.hardware.sequence_inst.ACQ.freq` and :meth:`~matipo.hardware.sequence_inst.ACQ.phase` commands.
    #. CIC Decimator: Decimates the quadrature data to the desired dwell time, set by the :meth:`~matipo.hardware.sequence_inst.ACQ.dwelltime` command. The antialiasing filter has a fast settling time of 6 samples but does not have a flat frequency response.
    #. FIR Compensation Filter: compensates for the response of the CIC filter to produce a flat frequency domain with a sharp rolloff. This filter has a settling time of approx. 20 samples.

The DDC may be disabled entirely by setting ``raw=True`` using the :meth:`~matipo.hardware.sequence_inst.ACQ.mode` command. In ``raw`` mode the output data will be the unprocessed 16 bit samples from the ADC.

.. note::
    In ``raw`` mode, all acquisitions must have a number of samples that is a multiple of 4.

The flat filter may be disabled by setting ``flatfilter=False`` using the :meth:`~matipo.hardware.sequence_inst.ACQ.mode` command. This may be desirable in applications that require fast settling time.

General Purpose Outputs
^^^^^^^^^^^^^^^^^^^^^^^

The GPO state can be controlled using the :meth:`~matipo.hardware.sequence_inst.GPO.gpo_set` and :meth:`~matipo.hardware.sequence_inst.GPO.gpo_clear` commands. Here is an example of using a GPO to control an active T/R switch while transmitting an RF pulse::

    TR_SWITCH_ENABLE = 0x00000001 # Port 3 pin 1 bitmask
    yield seq.gpo[0].set(TR_SWITCH_ENABLE) # set high, T/R switch transmit mode
    yield seq.tx[0].phase(0) # 0 degrees phase
    yield seq.tx[0].amp(0.5) # 50% amplitude
    yield seq.tx[0].enable() # enable RF output
    yield seq.wait(100e-6) # wait 100 microseconds
    yield seq.tx[0].disable() # disable RF output
    yield seq.gpo[0].clear(TR_SWITCH_ENABLE) # clear low, T/R switch receive mode

GPO Pin Mapping
"""""""""""""""

GPO Ports are indexed from left to right from the front, this means port 0 is next to the USB ports and port 3 next to the power port. There are 7 general purpose outputs on each port, using pin 1-6 and pin 8.

+-----+---------------------------------------------------+
|     |                     GPO Port                      |
+=====+============+============+============+============+
| Pin | 0          | 1          | 2          | 3          |
+-----+------------+------------+------------+------------+
| 1   | 0x00000001 | 0x00000100 | 0x00010000 | 0x01000000 |
+-----+------------+------------+------------+------------+
| 2   | 0x00000002 | 0x00000200 | 0x00020000 | 0x02000000 |
+-----+------------+------------+------------+------------+
| 3   | 0x00000004 | 0x00000400 | 0x00040000 | 0x04000000 |
+-----+------------+------------+------------+------------+
| 4   | 0x00000008 | 0x00000800 | 0x00080000 | 0x08000000 |
+-----+------------+------------+------------+------------+
| 5   | 0x00000010 | 0x00001000 | 0x00100000 | 0x10000000 |
+-----+------------+------------+------------+------------+
| 6   | 0x00000020 | 0x00002000 | 0x00200000 | 0x20000000 |
+-----+------------+------------+------------+------------+
| 8   | 0x00000040 | 0x00004000 | 0x00400000 | 0x40000000 |
+-----+------------+------------+------------+------------+

.. _error-codes:
Error Codes
^^^^^^^^^^^

If the pulse sequence execution driver detects an error, the pulse sequence will be terminated and the `Sequence.run()` method will raise an exception with an error code:

* -131: RFTX frequency out of bounds
* -132: RFTX amplitude out of bounds
* -141: ACQ frequency out of bounds
* -142: ACQ raw mode acquisition length not divisible by 4
* -144: ACQ dwell time out of bounds
* -145: ACQ number of samples out of bounds
* -146: ACQ data too large for DMA buffer
* -162: Sequencer wait time out of bounds (too small, must be >=100e-9 seconds)
* -171: SHIM channel out of range
* -1106: Input buffer ran out of pulse sequence commands before the end of the sequence. Usually may be remedied by optimising the sequence code, contact support for guidance.