Dear Peixiang, Dear Andrew.
Just a little heads-up.
Within a week or two, we should be able to release a new version of relax.
This update has focused on speed, recasting the data to higher
dimensional numpy arrays, and
moving the multiple dot operations out of the for loops.
So we have quite much better speed now. :-)
But we still have a bottleneck with numerical models, where doing the
matrix exponential via eigenvalue decomposition is slow!
Do any of you any experience with a faster method for doing matrix
exponential?
These initial results shows that if you are going to use the R1rho
2site model, you can expect:
-R1rho
100 single spins analysis:
DPL94: 23.525+/-0.409 -> 1.138+/-0.035, 20.676x
faster.
TP02: 20.191+/-0.375 -> 1.238+/-0.020, 16.308x
faster.
TAP03: 31.993+/-0.235 -> 1.956+/-0.025, 16.353x
faster.
MP05: 24.892+/-0.354 -> 1.431+/-0.014, 17.395x
faster.
NS R1rho 2-site: 245.961+/-2.637 -> 36.308+/-0.458, 6.774x
faster.
Cluster of 100 spins analysis:
DPL94: 23.872+/-0.505 -> 0.145+/-0.002, 164.180x
faster.
TP02: 20.445+/-0.411 -> 0.172+/-0.004, 118.662x
faster.
TAP03: 32.212+/-0.234 -> 0.294+/-0.002, 109.714x
faster.
MP05: 25.013+/-0.362 -> 0.188+/-0.003, 132.834x
faster.
NS R1rho 2-site: 246.024+/-3.724 -> 33.119+/-0.320, 7.428x
faster.
-CPMG
100 single spins analysis:
CR72: 2.615+/-0.018 -> 1.614+/-0.024, 1.621x
faster.
CR72 full: 2.678+/-0.020 -> 1.839+/-0.165, 1.456x
faster.
IT99: 1.837+/-0.030 -> 0.881+/-0.010, 2.086x
faster.
TSMFK01: 1.665+/-0.049 -> 0.742+/-0.007, 2.243x
faster.
B14: 5.851+/-0.133 -> 3.978+/-0.049, 1.471x
faster.
B14 full: 5.789+/-0.102 -> 4.065+/-0.059, 1.424x
faster.
NS CPMG 2-site expanded: 8.225+/-0.196 -> 4.140+/-0.062, 1.987x
faster.
NS CPMG 2-site 3D: 240.027+/-3.182 -> 45.056+/-0.584, 5.327x
faster.
NS CPMG 2-site 3D full: 240.910+/-4.882 -> 45.230+/-0.540, 5.326x
faster.
NS CPMG 2-site star: 186.480+/-2.299 -> 36.400+/-0.397, 5.123x
faster.
NS CPMG 2-site star full: 187.111+/-2.791 -> 36.745+/-0.689, 5.092x
faster.
Cluster of 100 spins analysis:
CR72: 2.610+/-0.035 -> 0.118+/-0.001, 22.138x
faster.
CR72 full: 2.674+/-0.021 -> 0.122+/-0.001, 21.882x
faster.
IT99: 0.018+/-0.000 -> 0.009+/-0.000,
2.044x faster. IT99: Cluster of only 1 spin analysis, since v. 3.2.2
had a bug with clustering analysis.:
TSMFK01: 1.691+/-0.091 -> 0.039+/-0.000, 43.369x
faster.
B14: 5.891+/-0.127 -> 0.523+/-0.004, 11.264x
faster.
B14 full: 5.818+/-0.127 -> 0.515+/-0.021, 11.295x
faster.
NS CPMG 2-site expanded: 8.167+/-0.159 -> 0.702+/-0.008, 11.638x
faster.
NS CPMG 2-site 3D: 238.717+/-3.025 -> 41.380+/-0.950, 5.769x
faster.
NS CPMG 2-site 3D full: 507.411+/-803.089 -> 41.852+/-1.317,
12.124x faster.
NS CPMG 2-site star: 187.004+/-1.935 -> 30.823+/-0.371, 6.067x
faster.
NS CPMG 2-site star full: 187.852+/-2.698 -> 30.882+/-0.671, 6.083x
faster.
The reason we cant tweak NS R1rho 2-site anymore, is that we are computing the
matrix exponential
##### Lib function for NS R1rho 2-site.
Before:
http://svn.gna.org/svn/relax/trunk/lib/dispersion/ns_r1rho_2site.py
Now: moving all essential computations to be calculated before.
http://svn.gna.org/svn/relax/branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py
def ns_r1rho_2site(M0=None, M0_T=None, .....
# Once off parameter conversions.
pB = 1.0 - pA
k_BA = pA * kex
k_AB = pB * kex
# Extract shape of experiment.
NE, NS, NM, NO = num_points.shape
# The matrix that contains all the contributions to the evolution,
i.e. relaxation, exchange and chemical shift evolution.
R_mat = rr1rho_3d_2site_rankN(R1=r1, r1rho_prime=r1rho_prime,
dw=dw, omega=omega, offset=offset, w1=spin_lock_fields, k_AB=k_AB,
k_BA=k_BA, relax_time=relax_time)
# This matrix is a propagator that will evolve the magnetization
with the matrix R.
Rexpo_mat = matrix_exponential_rank_NE_NS_NM_NO_ND_x_x(R_mat)
# Magnetization evolution.
Rexpo_M0_mat = einsum('...ij,...jk', Rexpo_mat, M0)
# Magnetization evolution, which include all dimensions.
MA_mat = einsum('...ij,...jk', M0_T, Rexpo_M0_mat)[:, :, :, :, :, 0, 0]
# Do back calculation.
back_calc[:] = -inv_relax_time * log(MA_mat)
#######
The profiling scripts, shows that
matrix_exponential_rank_NE_NS_NM_NO_ND_x_x(R_mat) is stealing the
time.
#####
http://svn.gna.org/svn/relax/branches/disp_spin_speed/lib/dispersion/matrix_exponential.py
matrix_exponential_rank_NE_NS_NM_NO_ND_x_x(A
# The eigenvalue decomposition.
W, V = eig(A)
# Calculate the exponential of all elements in the input array
# Add one axis, to allow for broadcasting multiplication
W_exp = exp(W).reshape(NE, NS, NM, NO, ND, Row, 1)
# Make a eye matrix, with Shape [NE][NS][NM][NO][ND][X][X]
eye_mat = tile(eye(Row)[newaxis, newaxis, newaxis, newaxis,
newaxis, ...], (NE, NS, NM, NO, ND, 1, 1) )
# Transform it to a diagonal matrix, with elements from vector down th
W_exp_diag = multiply(W_exp, eye_mat)
# Make dot products for higher dimension.
dot_V_W = einsum('...ij,...jk', V, W_exp_diag)
# Compute the (multiplicative) inverse of a matrix.
inv_V = inv(V)
# Calculate the exact exponential.
eA = einsum('...ij,...jk', dot_V_W, inv_V)
###
The numpy implementation of eig() is stealing up 86% of the time.
If you have any other way to do this efficient, I would be happy to hear it!
I have looked through this, do you have any experience with some of the
methods?
Moler, C. and Van Loan, C. (2003) Nineteen Dubious Ways to Compute
the Exponential of a Matrix, Twenty-Five Years Later. SIAM Review,
45, 3-49. (http://dx.doi.org/10.1137/S00361445024180 or
http://www.cs.cornell.edu/cv/researchpdf/19ways+.pdf).
Best
Troels
2014-06-11 10:36 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi Peixiang,
Actually, for comparison purposes for applying the 'NS R1rho 2-site'
model (http://wiki.nmr-relax.com/NS_R1rho_2-site) to variable-time
R1rho-type data, Art Palmer's MP05 model would be much better
(http://wiki.nmr-relax.com/MP05) than the DPL94 model
(http://wiki.nmr-relax.com/DPL94) as it is of much higher quality.
Andy Baldwin apparently has derived an even better analytic model,
especially when R20A and R20B are significantly different, see:
http://gna.org/support/?3155#comment0
and the discussions in the thread:
http://thread.gmane.org/gmane.science.nmr.relax.devel/5414/focus=5447
and:
http://thread.gmane.org/gmane.science.nmr.relax.devel/5410/focus=5433
This last thread is about the B14 model (Baldwin 2014,
http://wiki.nmr-relax.com/B14) implemented in relax by Troels Linnet,
but there are mentions of Andy's R1rho model. However the R1rho model
from Andy is not implemented in relax yet. Do you have much
experience with variable-time R1rho numeric models? Looking at the
code for where the relax_time variable comes from, it is not very
clear which relaxation time is being used:
http://www.nmr-relax.com/api/3.2/specific_analyses.relax_disp.data-module.html#loop_time
From the code itself:
http://www.nmr-relax.com/api/3.2/specific_analyses.relax_disp.data-pysrc.html#loop_time
it looks like this loop_time() function assumes fixed-time data and
hence only the first encountered time value for the given experiment,
magnetic field strength, offset, and dispersion point is used. So
your expertise will be very useful for resolving this variable-time
R1rho numeric model problem!
Note that there are a few improvements to the R1rho models that are
yet to be implemented in relax:
http://thread.gmane.org/gmane.science.nmr.relax.devel/5414/focus=5808
http://www.nmr-relax.com/manual/do_dispersion_features_yet_be_implemented.html
Cheers,
Edward
On 11 June 2014 10:07, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:
Hi Peixiang,
Please see below:
Congratulations about the new version of 3.2.2, I tried, it works well :)
Cheers. If you notice any other problems or strange behaviour, please
don't hesitate to submit a bug report
(https://gna.org/bugs/?func=additem&group=relax). Then that problem
will likely be fixed for the next relax version.
Still one question of using the different relaxation time periods.
My R1rho RD experiment has different relaxation time periods, I could
input all the peaks by the loop.
Then I fit with 'NS 2-site R1 model', they could also do the fitting and
give the results and also a nice fitting of the dispersion curve.
Still I did not figure out, which Trelax is it using in the NS model in
the case of different relaxation time periods.
Only the last relaxation time period? Then fit as fixed time experiment?
As this code was directly contributed by Paul Schanda and Dominique
Marion, and I'm guessing that their offices are not too far from yours
at the IBS, maybe you could ask them directly ;) Well, it was Paul
who organised that the code be contributed to relax. In reality the
original authors were Nikolai Skrynnikov and Martin Tollinger. The
API documentation is also a useful resource for answering such
questions (http://www.nmr-relax.com/api/3.2/). For this, see the
relax library documentation for that model:
http://www.nmr-relax.com/api/3.2/lib.dispersion.ns_r1rho_2site-module.html
This documentation describes the origin and history of the code. You
could even look at the source code for the direct implementation:
http://www.nmr-relax.com/api/3.2/lib.dispersion.ns_r1rho_2site-pysrc.html
Trelax is the 'relax_time' argument here. You can find all
implementation details in this API documentation. Which relaxation
time would you suggest as being correct? I'm actually no longer sure
which is being used. And I'm not sure if the original code or even
the numeric model itself was designed to handle variable time data.
Maybe I am the minority to use such time consuming experiments, so I
always have such strange questions ...
relax should still handle the situation. Do you know if there is a
special treatment for the numerical models for such data? Do you know
of a good citation? Maybe the 'NS R1rho 2-site' model
(http://wiki.nmr-relax.com/NS_R1rho_2-site) is not suitable for
variable time data, and a different - and importantly published -
solution is required. The analytic models do not use the relaxation
time value, so those are safe. Hence, as a check, you should see very
similar results from the 'DPL94' model
(http://wiki.nmr-relax.com/DPL94) and the 'NS R1rho 2-site' model. If
not, something is wrong.
If the 'NS R1rho 2-site' model is really only for fixed-time data,
then we should modify relax to raise a RelaxError when this model is
chosen for optimisation and the data is variable time. As not many
people optimise numeric models to variable-time data, your input into
this question would be very valuable. Cheers!
Maybe another annoying question for the fix time people:
Another question, does it necessary to check how mono-exponential about
their relaxation curve under certain rf-field? If not, how did they make
sure they can use the mono-exponential assumption to get R2eff by two
points?
From what I've seen and heard, some people do check, but the majority
just assume that the curves will be mono-exponential and publish the
fixed-time data and results. Such a check is probably much more
important for those collecting R1rho-type data rather than CPMG-type
data. Anyway, maybe you should ask people in front of their posters
at conferences to get a better overview of what the field does.
Regards,
Edward
Best,
Peixiang
On 05/19/2014 05:49 PM, Edward d'Auvergne wrote:
Hi Peixiang,
Welcome to the relax mailing lists! The relaxation dispersion
analysis implemented in relax is quite flexible, and the data you have
is supported. This is well documented in the relax manual which you
should have with your copy of relax (the docs/relax.pdf file). Have a
look at section 'The R2eff model' in the dispersion chapter of the
manual (http://www.nmr-relax.com/manual/R2eff_model.html),
specifically the 'Variable relaxation period experiments' subsection.
Unfortunately the sample scripts are all for the fixed time dispersion
experiments. However you could have a look at one of the scripts used
for the test suite in relax:
test_suite/system_tests/scripts/relax_disp/exp_fit.py
This script is run in the test suite to ensue that the data you have
will always be supported. There are many more scripts in that
directory which you might find interesting. The 'r1rho_on_res_m61.py'
script also involve an exponential fit with many different relaxation
time periods.
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-users mailing list
relax-users@xxxxxxx
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-users
Re: R1rho RD analysis