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.
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?
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.
Regards,
Edward
Re: CST branch