Re: CST branch -- March 22, 2011 - 09:19

Hi,

Please see below:

On 10 March 2011 18:55, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:

Hi,

Please see below:

On 10 March 2011 18:09, Pavel Kaderavek <pavel.kaderavek@xxxxxxxxx> wrote:
> Hi,
>
> Within this mail we need to discuss several related problems we faced when
> we started to think about the following patch. If it will be discussed
> separately the connection would be lost, therefore we decided to send it
> together. Of course we don't assume that everything will be put into the
> code in one huge patch, but it will be split into several parts.
>
> Originally, we wanted to discuss the selection of the equation in the
> function setup_equations in class Mf (maths_fns/mf.py). Currently, the
> equation for calculation physical constants and linear combination of
> spectral densities are selected just based on the type of relaxation rate
> (R1,R2,NOE).
>
>        for i in xrange(data.num_ri):
>             # The R1 equations.
>              if data.ri_labels[i] == 'R1':
>                 data.create_csa_func[i] = comp_r1_csa_const
>                 data.create_csa_grad[i] = comp_r1_csa_const
>                 data.create_csa_hess[i] = comp_r1_csa_const
>                 data.create_dip_jw_func[i] = comp_r1_dip_jw
>                 data.create_dip_jw_grad[i] = comp_r1_dip_jw
>                 data.create_dip_jw_hess[i] = comp_r1_dip_jw
>                 data.create_csa_jw_func[i] = comp_r1_csa_jw
>                 data.create_csa_jw_grad[i] = comp_r1_csa_jw
>                 data.create_csa_jw_hess[i] = comp_r1_csa_jw
>
>             # The R2 equations.
>             elif data.ri_labels[i] == 'R2':
>                 data.create_dip_func[i] = comp_r2_dip_const
>                 data.create_dip_grad[i] = comp_r2_dip_const
>                 data.create_dip_hess[i] = comp_r2_dip_const
>                 data.create_csa_func[i] = comp_r2_csa_const
>                 data.create_csa_grad[i] = comp_r2_csa_const
>                 data.create_csa_hess[i] = comp_r2_csa_const
>                 data.create_rex_func[i] = comp_rex_const_func
>                 data.create_rex_grad[i] = comp_rex_const_grad
>                 data.create_dip_jw_func[i] = comp_r2_dip_jw
>                 data.create_dip_jw_grad[i] = comp_r2_dip_jw
>                 data.create_dip_jw_hess[i] = comp_r2_dip_jw
>                 data.create_csa_jw_func[i] = comp_r2_csa_jw
>                 data.create_csa_jw_grad[i] = comp_r2_csa_jw
>                 data.create_csa_jw_hess[i] = comp_r2_csa_jw
>
>             # The NOE equations.
>             elif data.ri_labels[i] == 'NOE':
>                 data.create_dip_jw_func[i] = comp_sigma_noe_dip_jw
>                 data.create_dip_jw_grad[i] = comp_sigma_noe_dip_jw
>                 data.create_dip_jw_hess[i] = comp_sigma_noe_dip_jw
>                 data.create_ri[i] = calc_noe
>                 data.create_dri[i] = calc_dnoe
>                 data.create_d2ri[i] = calc_d2noe
>                 if data.noe_r1_table[i] == None:
>                     data.get_r1[i] = calc_r1
>                     data.get_dr1[i] = calc_dr1
>                     data.get_d2r1[i] = calc_d2r1
>                 else:
>                     data.get_r1[i] = extract_r1
>                     data.get_dr1[i] = extract_dr1
>                     data.get_d2r1[i] = extract_d2r1
>
> It is necessary to select the equation based on both type of relaxation rate
> and also interaction type in the new code design. Therefore we suggest to
> change the code in following way:
>
>         for i in xrange(data.num_ri):
>             # The R1 equations.
>             if data.ri_labels[i] == 'R1':
>                 if data.interactions == "dip"
>                    data.create_const_func[i] = comp_dip_const_func
>                    data.create_const_grad[i] = comp_dip_const_grad
>                    data.create_const_hess[i] = comp_dip_const_hess
>                    data.create_jw_func[i] = comp_r1_dip_jw
>                    data.create_jw_grad[i] = comp_r1_dip_jw
>                    data.create_jw_hess[i] = comp_r1_dip_jw
>                 elif data.interactions == "CSA"
>                    data.create_const_func[i] = comp_csa_const_func
>                    data.create_const_grad[i] = comp_csa_const_grad
>                    data.create_const_hess[i] = comp_csa_const_hess
>                    data.create_jw_func[i] = comp_r1_csa_jw
>                    data.create_jw_grad[i] = comp_r1_csa_jw
>                    data.create_jw_hess[i] = comp_r1_csa_jw
>                 elif data.interactions == "Rex"
>                    data.create_const_func[i] = comp_rex_const_func
>                    data.create_const_grad[i] = comp_rex_const_grad
>                    data.create_const_hess[i] = comp_rex_const_hess
>                 elif data.interactions == "cross-CSA-CSA"
>                    data.create_const_func[i] = comp_cross_csa_csa_const_func
>                    data.create_const_grad[i] = comp_cross_csa_csa_const_grad
>                    data.create_const_hess[i] = comp_cross_csa_csa_const_hess
>                    data.create_jw_func[i] = comp_r1_cross_csa_csa_jw
>                    data.create_jw_grad[i] = comp_r1_cross_csa_csa_jw
>                    data.create_jw_hess[i] = comp_r1_cross_csa_csa_jw
>
>             # The R2 equations.
>             elif data.ri_labels[i] == 'R2':
>                 if data.interactions == "dip"
>                    data.create_const_func[i] = comp_dip_const_func
>                    data.create_const_grad[i] = comp_dip_const_grad
>                    data.create_const_hess[i] = comp_dip_const_hess
>                    data.create_jw_func[i] = comp_r2_dip_jw
>                    data.create_jw_grad[i] = comp_r2_dip_jw
>                    data.create_jw_hess[i] = comp_r2_dip_jw
>                 elif data.interactions == "CSA"
>                    data.create_const_func[i] = comp_csa_const_func
>                    data.create_const_grad[i] = comp_csa_const_grad
>                    data.create_const_hess[i] = comp_csa_const_hess
>                    data.create_jw_func[i] = comp_r2_csa_jw
>                    data.create_jw_grad[i] = comp_r2_csa_jw
>                    data.create_jw_hess[i] = comp_r2_csa_jw
>                 elif data.interactions == "Rex"
>                    data.create_const_func[i] = comp_rex_const_func
>                    data.create_const_grad[i] = comp_rex_const_grad
>                    data.create_const_hess[i] = comp_rex_const_hess
>                 elif data.interactions == "cross-CSA-CSA"
>                    data.create_const_func[i] = comp_cross_csa_csa_const_func
>                    data.create_const_grad[i] = comp_cross_csa_csa_const_grad
>                    data.create_const_hess[i] = comp_cross_csa_csa_const_hess
>                    data.create_jw_func[i] = comp_r2_cross_csa_csa_jw
>                    data.create_jw_grad[i] = comp_r2_cross_csa_csa_jw
>                    data.create_jw_hess[i] = comp_r2_cross_csa_csa_jw
>
>             # The NOE equations.
>             elif data.ri_labels[i] == 'NOE':
>                 if data.interactions == "dip"
>                    data.create_const_func[i] = comp_dip_const_func
>                    data.create_const_grad[i] = comp_dip_const_grad
>                    data.create_const_hess[i] = comp_dip_const_hess
>                    data.create_jw_func[i] = comp_sigma_noe_dip_jw
>                    data.create_jw_grad[i] = comp_sigma_noe_dip_jw
>                    data.create_jw_hess[i] = comp_sigma_noe_dip_jw
>
>                 data.create_ri[i] = calc_noe
>                 data.create_dri[i] = calc_dnoe
>                 data.create_d2ri[i] = calc_d2noe
>                 if data.noe_r1_table[i] == None:
>                     data.get_r1[i] = calc_r1
>                     data.get_dr1[i] = calc_dr1
>                     data.get_d2r1[i] = calc_d2r1
>                 else:
>                     data.get_r1[i] = extract_r1
>                     data.get_dr1[i] = extract_dr1
>                     data.get_d2r1[i] = extract_d2r1

This seems to be the correct way to go.  This setting up function is
quite complex, but it does save incredible amounts of computation
time.  So complexity here is not an issue.  With proper system and
unit testing, any problems should be quick to locate.  Speaking of a
system test, one should be designed for this new cst branch so that we
know when the code is fully functional.

Alright, so we will prepare next patch with the above changes.

> At this point we would like to address a related question. Currently the
> calculation of physical constant is done in a several steps. First, the
> physical constant is calculated and the value is stored in the
> data.dip_const_func or data.csa_const_func (grad, hess). Then, when the
> relaxation rates are calculated, the physical constant is recalculated by
> the function create_dip_func or create_csa_func (grad, hess) (method
> setup_equations in class Mf, maths_fns/mf.py).
>
>             comp_dip_const_func(data, data.bond_length)
>             comp_csa_const_func(data, data.csa)
>             for i in xrange(data.num_ri):
>                 data.dip_comps_func[i] = data.dip_const_func
>                 if data.create_dip_func[i]:
>                     data.dip_comps_func[i] =
> data.create_dip_func[i](data.dip_const_func)
>                 if data.create_csa_func[i]:
>                     data.csa_comps_func[i] =
> data.create_csa_func[i](data.csa_const_func[data.remap_table[i]])
>
> There is one exception, the dipolar physical constant is not recalculated in
> the case of calculation R1 relaxation rate, because the function
> create_dip_func does not exist in this case. We do not see a reason for such
> a recalculation.

The reason is because of the m10 to m39 models built into relax.  I
have made it possible to optimise the bond length and CSA information.
 However these models are not stable with the current relaxation data.
 I do plan on working with these in the future though, so it would be
useful to keep them.  Note that for models m0 to m9, the
data.create_dip_func[i] and data.create_csa_func[i] function pointers
are set to None.  Therefore for normal model-free analysis the
constant is not recalculated.

> It seems better to us to just change the coefficient in the
> functions comp_r1_dip_jw, comp_r2_dip_jw, comp_r1_csa_jw, comp_r2_csa_jw
> (maths_fns/ri_comps.py). I guess, that this design was dedicated to avoid
> multiple calculation of the same interaction constant for each measured
> relaxation rate. We would suggest to reach the same effect by this
> construction:
>
>             for i in xrange(data.num_ri):
>                 if data.const_func[0]:
>                     data.const_func[i] = data.const_func[0]
>                 else
>                     data.create_const_func(data)

For models m10 to m39, this construct will not work.  The constants
are already pre-calculated for models m0 to m9 so this is not needed.

> Note, the comp_dip_const_func and comp_csa_const_func should be change so
> that, it is possible to call them just with the argument data
> (maths_fns/ri_comps.py). Instead of:
>
>
> def comp_dip_const_func(data, bond_length):
>     """Calculate the dipolar constant.
>
>    ...
>
>     if bond_length == 0.0:
>         data.dip_const_func = 1e99
>     else:
>         data.dip_const_func = 0.25 * data.dip_const_fixed * bond_length**-6
>
>
> It should look like:
>
> def comp_dip_const_func(data):
>     """Calculate the dipolar constant.
>
>    ...
>
>     if data.internuclei_distance == 0.0:
>         data.const_func = 1e99
>     else:
>         data.const_func = 0.25 * data.dip_const_fixed *
> data.internuclei_distance**-6

The bond_length arg was designed so that the bond length could either
come from a fixed value supplied by the user or from the parameter
vector when bond lengths or CSA values are optimised.  This behaviour
might have to be preserved.

> data.dip_const_func were renamed to more general data.const_func and instead
> of bond_length the function directly takes  the internuclei distance for the
> current dipole-dipole interaction. The change of data.dip_const_func to
> data.const_func later simplify the code design in the maths_fns/ri_prime.py
> . It will be reduced just to a multiplication of constant and the linear
> combination of spectral density functions.

For models m10 to m39, I'm not sure if this design would work.  Could
we redesign this in another way in which these complex models are
still functional?

Then, we would suggest to call the function  comp_dip_const_func, comp_csa_const_func ... with full set of possible parameters, i.e.
comp_dip_const_func(data,internuclei_distance,csa1,csa2,rex)
comp_csa_const_func(data,internuclei_distance,csa1,csa2,rex)
...
If we would call the function with just physically relevant arguments then we would have to use yet another condition to decide the type of the interaction in the loop where individual interaction contributions are calculated. The physical quantities irrelevent for the given interaction (for example csa1,csa2,rex for the dipole-dipole interaction) are None from the initialization and will not be used by the function anyway.

> Moreover, there is an unanswered question about the NOE and the additional
> dipolar interaction. I am not sure if the suggested design is physically
> correct, rather not. During the NOE experiment, the protons are saturated in
> order to reach the steady state. Then a complete set cross relaxation rates
> between all interacting spin pairs should be taken into the account, not
> only between the spin of interest and all other interacting nuclei. On the
> other hand this is probably beyond the aim of the program relax. What do you
> think about that?

This is getting quite complex as you would need to take the
cross-correlated relaxation rates between the different relaxation
interactions into account, as well as the motion of all spins if they
are not directly bonded.  Is this needed for the current work?  Of
course anything is accepted into relax, especially if you would like
to probe this area (with a paper in mind), but it has to play nicely
with the rest of relax and not be a burden on the relax developers to
maintain in the future (as well as not make the current number
crunching code slower than it already is).  The code would therefore
need to be designed in public.  So if you would like to tackle such a
task, I would first recommend finishing off the cst branch, and then
make a new branch for this work.

Then, I will design the code as I suggested. So, the sigma_noe will be calculated separately for all dipole-dipole interactions with the central spin (assuming they are isolated spin pairs). The total sigma_noe will be calculated as a sum of all individual sigma_noe contributions.

Regards,

Edward

Best regards,
Pavel