The following is a tutorial for adding new relaxation dispersion
models for either CPMG-type or R1rho-type experiments to relax. This
includes both the models based on the analytical, closed-form
expressions as well as the models involving numerical integration of
the Bloch-McConnell equations.
The tutorial will follow the example of the addition of the 'M61'
model already present within relax, pointing to the relevant commits
for reference. To see the commit message and the code changes in
colour, click on the links found within these commit messages. This
specific case is the Meiboom 1961 analytic model for 2-site fast
exchange equation for R1rho-type experiments.
1) Adding the model to the list.
Reference commit: http://article.gmane.org/gmane.science.nmr.relax.scm/17611
Firstly the model should be added to the lists of the
specific_analyses.relax_disp.variables module. The model name is
stored in a special variable which will be used throughout relax.
2) The relax_disp.select_model user function front end.
Reference commit: http://article.gmane.org/gmane.science.nmr.relax.scm/17612
The next step is to add the model, its description, the equations for
the analytic models, and all references to the relax_disp.select_model
user function front end. When the relaxation dispersion chapter of
the relax manual is created (this will be the
docs/latex/relax_disp.tex file), then the same description should be
added there as well.
3) The relax library.
Reference commit: http://article.gmane.org/gmane.science.nmr.relax.scm/17615
Now the dispersion function needs to be added to the relax library (in
the lib.relax_disp package). This should be designed as a simple
Python function which takes the dispersion parameters and experimental
variables, and calculates the R2eff/R1rho values. The module can
contain auxiliary functions for the calculation. Some auxiliary
functions, if not specific to relaxation dispersion, may be better
placed in other locations within the relax library.
The relaxation dispersion functions in the library currently take as
an argument a data structure for the back-calculated R2eff/R1rho
values and populate this structure. This design is not essential if
the target function, described in the next point, handles the library
function appropriately. Just look at the files in lib/dispersion to
get an idea of the design used.
The dispersion code in the relax library must be robust. This
involves identifying parameter values or combinations which would
cause failures in the mathematical operations (numerical issues not
present in the mathematics must be considered). Note that parameter
values of 0 are common within a grid search. It should be decided if
the R2eff/R1rho value should be set to zero, to another value, or to
something large (e.g. 1e100). For example:
Divisions - always catch zeros in the denominator with if statements,
even if you believe that this will never be encountered.
Square roots - make sure that the value inside is always > 0.
Trigonometric functions - these should be tested for where they are
not defined or where the software implementation can no longer handle
certain values. For example try cosh(1000) in Python.
In the reference example, the M61 model code was copied from the LM63
module and modified appropriately.
4) The target function.
Reference commits:
http://article.gmane.org/gmane.science.nmr.relax.scm/17616
http://article.gmane.org/gmane.science.nmr.relax.scm/17660
http://article.gmane.org/gmane.science.nmr.relax.scm/17661
The target function is used in optimisation and is a class method
which takes as a single argument the parameter vector. This list is
changed by the minimisation algorithm during optimisation. The target
function should then return a single floating point number - the
chi-squared value.
Again in this example, the code for the M61 is copied from the LM63
model and then modified.
5) Adding support for the parameters.
Reference commit: http://article.gmane.org/gmane.science.nmr.relax.scm/17573
This is needed to enable the model. This example is for the CR72
model implementation as the parameters required for the M61 model
match those of the preexisting LM63 model.
6) The relax_disp.select_model back end.
Reference commit: http://article.gmane.org/gmane.science.nmr.relax.scm/17622
Now the back end of the relax_disp.select_model user function for the
model can be added. This involved identifying the model and
constructing the parameter list.
7) The test suite.
Reference commits:
http://article.gmane.org/gmane.science.nmr.relax.scm/17647
http://article.gmane.org/gmane.science.nmr.relax.scm/17648
http://article.gmane.org/gmane.science.nmr.relax.scm/17649
http://article.gmane.org/gmane.science.nmr.relax.scm/17662
http://article.gmane.org/gmane.science.nmr.relax.scm/17663
This step is normally performed as step number 1. This is the most
important part that makes sure that the code not only works now, but
will continue working for the entire lifetime of the relax project.
The idea is that synthetic data (here for example as Sparky peak
lists) is created for the model and added to the test suite directory
test_suite/shared_data/dispersion/. This is then used in a system
test to check that the code in relax can reproduce the data. It is
very important that the code added to the relax library is not used to
create the synthetic data!
8) Comparing to other software.
It can happen that a bug present in the lib.dispersion package code is
also replicated in the synthetic data. This is not uncommon.
Therefore it is very useful to use other software with the test data
from step 7 to see if the original parameters can be found. A good
example can be seen in the test_suite/shared_data/dispersion/Hansen
which contains Dr. Flemming Hansen's CPMG data (see the README file)
and the results from different programs including NESSY, relax,
CPMGFit, and ShereKhan. The comparison is in the file
'software_comparison'.
Once the relax code is able to find identical or better results than
the dispersion softwares, then the values found in the test suite
optimisation can be locked in. The assertEqual() and
assertAlmostEqual() methods can be used to only allow the test to pass
when the correct values are found.
9) Debugging.
This step should not require an explanation. It goes hand-in-hand
with steps 7) and 8).
Tutorial for adding relaxation dispersion models to relax.