A not so Short Tutorial on DAMARIS

Python Front End for Experimenters

DAMARIS is the software, first written and conceived by Dr. Achim Gädke, used to control the NMR spectrometers in our lab. There are two separate parts: the back end, which controls the hardware and executes the pulse program, and a front end where the pulse program is defined and results are displayed. The programming language used for the scripts in the front end is Python, because it is easy to learn and the programs are easy to maintain. A introduction in Python is out of the scope for this thesis, however, because Python is widely used outside DAMARIS, there are a wealth of excellent tutorials available.

Basics

The basic idea in DAMARIS is the separation of an experiment in three steps:

  1. Experiment description where the experiment sequence is described
  2. Back end executes the sequences obtained from the experiment description
  3. Results from the back end are processed in the result script

This introduction into DAMARIS will begin with a NMR “Hello World” program, i.e. recording a free induction decay. More elaborated techniques like phase cycling and PFG will be discussed later.

Frontend User Interface

The front end is composed of five tabs (Figure 1.1): Experiment, Result, Display, Log and Configuration.

The first step is the configuration of DAMARIS. This is done with in th configuration tab which offers a convenient way to setup DAMARIS. Following is a short description of the options:

  1. Toggle the start of the back end
  2. Path to the back end executable, here PFGcore.exe
  3. Path to the spool directory, where the job files will be written or picked up
  4. Delete the job files after execution
  5. Filename of the logfile
  6. Number of logfiles to keep
  7. Toggle the start of the experiment script
  8. Toggle the start of the result script
  9. Toggle the deletion of results after they have been processed by the result script (see 1.8)
  10. Filename of the data pool file
  11. Write interval after which the data pool is rewritten
  12. Compression settings (see page 20)
  13. Authors, License and version numbers of used components

Experiment Script

Let us start with the most simple NMR experiment: recording a free induction decay, or short FID. This means we apply a RF pulse and record the decay of the sample magnetization. The corresponding DAMARIS experiment script would look like this:

def fid_experiment():
	    # create a sequence
	    e=Experiment()
	    e.set_frequency(frequency=300.01e6, phase=0)# set the frequency
	    e.ttl_pulse(length=5e-6, channel=1)    # gate for the RF pulse
	    e.ttl_pulse(length=2e-6, channel=1+2)    # RF pulse
	    e.wait(10e-6) # dead time
	    # record 1024 samples at sampling-rate 2 MHz
	    # with +/- 2V range
	    e.record(samples=1024, frequency=2e6, sensitivity=2)
	    return e

	def experiment():
	    yield fid_experiment()

On the first line, the basic Experiment module is imported which loads the methods and commands to control an experiment. Next, a function fid_experiment is defined which will contain everything necessary to run a single scan. The first line inside the function creates an object e in which the states will be written by subsequent application of methods to this experiment object. This sequence of states is then passed to the caller. When the “Execute”  button is pressed, the front end will execute a function called experiment. Thus, to run fid_experiment, the experiment function has to contain a call to the fid_experiment.

Behind the scenes a state sequence from fid_experiment, an XML2 file, will be written into the spool directory and the back end is then executed. The back end simply reads these files, translates the contents to pulse card command code and loads the generated code into the memory of the pulse card. Finally, the pulse program is executed. After succesful execution the results are written by the back end into the spool directory, where they are provided to the result script. Finally, in the result script, the function result is then executed for each obtained result.

Result Script

The function result in the result tab defines what should be done with the results. Here is a  very short example how to display the recorded FID from our experiment:

def result():
	    # loop over the incoming results
	    for timesignal in results:
	        # provide the timesignal to the display
	        data["Timesignal"]=timesignal

After pressing the “Execute” () button to start the experiment, the Display tab will open. With the drop down menu  under the graph window  (to the right of “Source”) one can choose  the data which should be displayed. This list is built dynamically by the data dictionary in the result script, where the key denotes the name to be shown and the value contains the data. In this example we overwrite the data with the key "Timesignal" on every scan with the new result. There are now two possibilities to export the figure:

  1. Saving the plot as a figure (PDF, EPS, PNG, JPEG and more) by pressing the button   in the toolbar below the plot area
  2. Saving the data of the plot as CSV file (comma separated values) by pressing the button  on the right side of the toolbar under the plot

Accumulation

Sometimes the signal intensity is low compared to the noise, and can therefore not be seen. A measure for the strength of the recorded signal is the signal-to-noise ratio (SNR) which is defined as the ratio of signal versus the standard deviation of the noise, i.e. noise strength. The SNR can be improved by measuring the signal many times and averaging the recorded signals, also called accumulation. In DAMARIS this can be achieved by creating a loop in the experiment script and accumulating the results in the result script.

def fid_experiment():
	    e=Experiment.Experiment()
	    e.set_frequency(frequency=300.01e6, phase=0)
	    e.ttl_pulse(length=5e-6, channel=1)  # gate for the RF pulse
	    e.ttl_pulse(length=2e-6, channel=1+2)  # RF pulse
	    e.wait(10e-6) # dead-time
	    # record 1024 samples at sampling-rate 2 MHz
	    # with 2V sensitivity
	    e.record(samples=1024, frequency=2e6, sensitivity=2)
	    return e
	
	def experiment(): 
            accumulations=8
            for run in xrange(accumulations):
                yield fid_experiment()

Our experiment script is hardly changed, we only added a  Python loop which will execute our experiment 8 times. The result script needs also only a minor adjustment:

	def result():
           # create accumulation object 
    accu = Accumulation.Accumulation()
     # loop over the incoming results 
    for timesignal in results: 
        # plot the timesignal 
        data["Timesignal"] = timesignal 
        # add timesignal to the object 
        accu += timesignal 
        # provide the data to display tab 
        data["Accumulation"] = accu

 

Phase Cycling

Artefacts due to channel imbalance in the receiver, inhomogeneity in the RF field and unwanted signals, for example primary echos, can be suppressed or cancelled out by cycling the phases of the receiver and the RF pulses and then accumulating the signal according to a phase cycling scheme. A very common scheme is the CYCLOPS phase cycling scheme, which cancels receiver channel imbalance and non othogonality of the real and imagninary part of the signal.


def fid_experiment(run): 
    e=Experiment() 
    # set a description, this one is used later in the result script 
    e.set_description("run",run)
    # change phase, CYCLOPS phase cycling
    # run%4 is run modulo(4), and thus goes like 0,1,2,3,0,1,2 ... 
    e.set_frequency(frequency=300.01e6, phase=[0,90,180,270][run%4]) 
    e.ttl_pulse(length=5e-6, channel=1)  # gate for the RF pulse 
    e.ttl_pulse(length=2e-6, channel=1+2)  # RF pulse 
    e.wait(10e-6)  # dead-time
    e.set_phase(0)
    # record 1024 samples at sampling-rate 2 MHz 
    # with a +/- 2V range 
    e.record(samples=1024, frequency=2e6, sensitivity=2) 
    return e def experiment(): 
    # should be multiple of 4 
    accumulations=8 
    for run in xrange(accumulations): 
        yield fid_experiment(run)

New in this script is the set_description method and a parameter run that is passed to the experiment. The set_description method is used to set descriptions and values important to the experimenter or to pass values to the result script. These descriptions are written into the XML job file and stored by the back end into the result file. In this way it is possible to pass values to the result script so the latter can act upon these. In this case the result is either added or subtracted to the accumulation variable.

 

Variable Parameters

Of course we want to change some parameters in order to gain insight of the dynamic or structure of the sample. After all this is what we want to do with NMR. There are a multidude of experiment where one or more parameters are varied in a certain way in order to elucidate some physical information, sometimes multiple parameters are variied which can result in  multi dimensional data.

As an simple yet typical exmaple, our experiment script is modified to measure the longitudinal relaxation time T1 via the so called inversion recovery (IR) experiment.In the IR experiment the 90°-pulse is preceded by a  180°-pulse both separated by a variable time here called tau. The FID after the second pulse is recorded and the signal height is analyzed. In order to get a good signal, accumulation and phase cycling are also implemented. In the result script, the magnetization will be obtained and, together with the corresponding variable τ, stored in an ASCII file for further analysis of T1, for exmaple fitting, etc.

Following is the experiment script and the result script which facilitate all of the above mentioned variation of parameters.


def fid_experiment(pi, tau, run, accus): 
    e=Experiment() 
    e.wait(9) # repetition time (in s) 
    # set the descriptions, they are used later in the result script 
    e.set_description("run",run) 
    e.set_description("tau",tau) 
    e.set_description("no_accus", accus) 
    # run%4 is run modulo(4), and thus goes like 0,1,2,3,0,1,2 ... 
    e.set_frequency(frequency=300.01e6, phase=[0,180][run%2]) 
    e.ttl_pulse(length=5e-6, channel=1)  # gate for the RF pulse 
    e.ttl_pulse(length=pi, channel=1+2)  # RF pulse 
    e.wait(tau) 
    e.set_phase(phase=[0,0,180,180,90,90,270,270][run%8]) 
    e.ttl_pulse(length=5e-6, channel=1)  # gate for the RF pulse 
    e.ttl_pulse(length=pi/2, channel=1+2)  # RF pulse 
    e.wait(10e-6)   # dead-time 
    # change phase, CYCLOPS phase cycling 
    e.set_phase([0,0,180,180,90,90,270,270][run%8]) 
    # record 1024 samples at sampling-rate 2 MHz 
    # with a +/- 2V range 
    e.record(samples=1024, frequency=2e6, sensitivity=2) 
    return e

def experiment(): 
    # should be multiple of 8 
    accumulations=8 
    # here we loop over the tau values starting from 1ms to 10s 
    # in a logarithmic 
    for t in log_range(start=1e-3, stop=10, step_no=20): 
        # this inner loop does the accumulation
        for r in xrange(accumulations): 
            yield fid_experiment(pi=4e-6, tau=t, run=r, accus=accumualations)

And here follows the corresponding result script:


def result(): 
    # create accumulation object 
    accu = Accumulation() 
    # loop over the incoming results 
    for timesignal in results: 
        # read in variables 
        tau = float(timesignal.get_description("tau")) 
        run = int(timesignal.get_description("run")) 
        accus = int(timesignal.get_description("no_accus"))-1 
        # provide the single scan to the display tab 
        data["Timesignal"] = timesignal 
        accu += timesignal 
        # provide the accumulated signal to display tab 
        data["Accumulation"] = accu 
        if run == accus: 
            # provide the accumulated signal to display tab
            # in this case all accumulated signals will be 
            # available in the data dictionary and thus 
            # saved when the experiment is finnished 
            data["Accumulation %i"%run] = accu 
            # open a file in append mode, 
            # file will be created if not existing 
            afile.open(’t1.dat’,’a’) 
            # get the maximum on the y channel 
            amplitude = accu.y[0].max() 
            # write tau and amplitude seperated by a 
<tab> and             # ending with a <newline> into the file 
            afile.write("%e\t%e\n"%(tau, amplitude)) 
            # close the file 
            afile.close() 
            # create a new accumulation object 
            accu = Accumulation()

Data Handling

For further data analysis one has several built-in possibilities in DAMARIS. The results can be exported as ASCII files or as HDF5 file (since 0.13svn also as SIMPSON file). Writing an ASCII file is usable for smaller data-sets and a smaller number of experiments. The big advantage is the broad usability: almost every program can import this data.

timesignal.write_to_csv(’filename_to_store.dat’)

For more elaborate experiments it can be advantageous to store data in a HDF5 file. This is the Hierarchical Data Format developed by the NCSA (National Center for Supercomputer Applications) for storing complex tables and incredible amounts of data in an effective way. It supports grouping of data, annotation and compression. A very convenient way to access HDF files in Python is the pyTables module, which is also used by DAMARIS internally. A platform independent way to examine the results is the HDFView Java program provided by NCSA.

After an experiment in DAMARIS has been finished, the results in the data dictionary as well as the experiment and result scripts themselves are written automatically to a HDF file, the data pool file. In the configuration tab  you can set a filename, compression ratio, compression library and write interval. The write interval is the time after which the data pool file should be periodically rewritten and is useful in case of hardware failure. The results from the ADC card are stored as 64-bit floating point values in tables which are compressed/decompressed transparently/on-the-fly. In DAMARIS, three compression libraries are available from the underlying pyTables mpdule: zlib, lzo and bz2. For a comparison of these, the pyTables homepage offers a great wealth of information. HDFView does not support viewing of bzip2 or lzo compressed data. In case the data has been compressed with lzo, the ptrepack utility from pyTables can be used to recompress the data with zlib.

The speed-up in post-processing and analysis is already 24 fold in small data-sets with 32768 points per channel and 25 files on a Athlon XP 2600+ CPU with 512 MB RAM running Debian Sarge Linux. Thus, writing HDF5 with compression enabled is highly recommended if your analysis software can read HDF5 files.

In the following code we will create a HDF5 file from scratch just to show what can be done! Experiment and result script are not saved in the resulting file. The advantage is that one can create arbitrary complex structured files at the cost of having a more complicated result script.


import tables 
def fid_result(): 
    # open file in "w"rite mode to generate a new, empty hdf5 file 
    # this is not absolutely necessary if the file does not already exists 
    nmr_file = tables.openFile(’filename_to_store.h5’, ’w’) 
    # open file in "a"ppend mode 
    nmr_file = tables.openFile(’filename_to_store.h5’, ’a’) 
    # loop over the results 
    for timesignal in results: 
        # plot the timesignal 
        data["Timesignal"]=timesignal 
        # write data 
        timesignal.write_to_hdf(hdffile=nmr_file, 
                            # save it under the root group of the file 
                            where="/", 
                            # name of the group 
                            name="accu", 
                            title="accudata", 
                            # enable lzo compression 
                            level="5", 
                            complib="zlib") 
    # flush the data to the disk 
    nmr_file.flush() 
    # close the file 
    nmr_file.close()
def result(): 
    return fid_result()

A very convenient and much easier method to save data in HDF5, or more precisely, in the data pool file is the MeasurementResult class.


def fid(): 
    # create a measurement object 
    measurement=MeasurementResult("Magnetization(tau)") 
    
for timesignal in results: 
    # extract "tau" from the description dictionary 
        tau = timesignal.get_description(’tau’) 
        ... 
        # get the magnetization from the timesignal 
        magnetization = sqrt(timesignal.y[0][80:100].mean()**2 
                        + timesignal.y[1][80:100].mean()**2) 
        measurement[tau] += magnetization 
        # provide the magnetization curve 
        data["Magnetization(tau)"] = measurement

Essentially, a dictionary “measurement” is created with tau as the key. For each new key the current magnetization is added. This dictionary is then added to the data pool. The advantage over the other code is that the magnetization vs. τ curve can be directly displayed in order to gain an overview of the running experiment, because the “measurement” dictionary is added to the data pool. In this example, the value also stores statistics from the accumulation (which does not make sense due to the phase cycling). The data dictionary is stored automatically as described before, including the result and experiment scripts, information about the used back end and a timeline of the experiment.

Good Practices

  • Keep the experiment function simple. Everything necessary for a scan should be contained in a function, the parameters should stay in the experiment function. This makes it possible to loop over several parameters easier and, in case of complex experiments, all variables are changed in one place: the call to the function.
  • Display only if necessary. If the repetition time of subsequent scans is short (ms) it can slow down the result script notably. The result script will take care of this by not displaying every single shot, but displaying nothing would further speed up the processing of the results if you are in a hurry.
  • Save disk space. If the options “delete results” and “delete jobs” are switched off, the spool directory can get clobbered by several thousand files and in the worst case, one can run out of inodes. If it is not possible to delete the files in the directory anymore (rm gives the error “too much arguments”) the following command can help deleting the files anyway: find ./spool -name job\* -exec rm {} \; which will delete all files starting with "job" in the spool directory.
  • It is a good idea to put a synchronize() statement in the experiment script, called for example after every 100th accumulation. This will write job files only untilsynchronize() is requested. Then, the front end will wait with writing new job files until the result script has processed all current results. In big experiments, several ten thousands job of files could be written without this, which can lead to problems with the file system and performance. if accu%100 == 0: synchronize()
  • Keep an eye on the logs. In case of errors it is highly recommended to check the log file in the log path (as defined in the configuration tab) as well as the “Log” tab for error messages.

Tuning and Development

When new result scripts are developed, it ican be sometimes an advantage if the actual experiment has not to be repeated for every minor change or tweak in the result script. This can become a very tedious task for example when the sample has a long T1 relaxation, like water T13s[4]. In order to use this technique, the experiment has to be run once with the option “delete results after processing” turned off. This leads to the result files in the spool directory not being consumed by the result script after processing them. Starting the experiment again with the option “run back end” switched off, the front end will read the results again like if it would be a new experiment, hence one can test new result scripts quite conveniently. Developing experiment scripts can be done in a similar manner: the dummy back end (to be set in the configuration tab) can be used to test new experiment scripts without the need to have a spectrometer at hand.

Data Processing Outside DAMARIS

For post-processing of the data, several options are available, depending on the format of the data. Most programs can handle ASCII files (gnuplot, Origin, fitsh, etc.), so this is the more flexible format. HDF can be read by Python with the help of pyTables, and h5py, in Matlab, Octave, R, Mathematica and many more.