Author: tlinnet
Date: Wed Jun 25 02:14:32 2014
New Revision: 24292
URL: http://svn.gna.org/viewcvs/relax?rev=24292&view=rev
Log:
Inserted an extremely interesting development in NS R1rho 2site.
If one do a transpose of M0, one can calculate all the matrix evolutions in
the start via numpy einsum.
Since M0 is in higher a dimensions, one should not do a numpy transpose,
but swap/roll the outer M0 6x1 axis.
Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion
models for Clustered analysis.
Modified:
branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py
Modified: branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py?rev=24292&r1=24291&r2=24292&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py (original)
+++ branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py Wed Jun 25
02:14:32 2014
@@ -50,7 +50,7 @@
"""
# Python module imports.
-from numpy import array, cos, dot, einsum, float64, log, multiply, sin, sum
+from numpy import array, cos, dot, einsum, float64, log, multiply,
rollaxis, sin, sum
# relax module imports.
from lib.float import isNaN
@@ -238,31 +238,63 @@
# Magnetization evolution.
Rexpo_M0_mat = einsum('...ij,...jk', Rexpo_mat, M0)
- # Loop over spins.
- for si in range(NS):
- # Loop over the spectrometer frequencies.
- for mi in range(NM):
- # Loop over offsets:
- for oi in range(NO):
- # Extract number of points.
- num_points_i = num_points[0, si, mi, oi]
-
- # Loop over the time points, back calculating the R2eff
values.
- for j in range(num_points_i):
+ # Transpose M0, to prepare for dot operation. Roll the last axis one
back.
+ M0_T = rollaxis(M0, 6, 5)
+
+ if NS != 1:
+ # Magnetization evolution, which include all dimensions.
+ MA_mat = einsum('...ij,...jk', M0_T, Rexpo_M0_mat)
+
+ # Loop over the spectrometer frequencies.
+ for mi in range(NM):
+ # Loop over offsets:
+ for oi in range(NO):
+ # Extract number of points.
+ num_points_i = num_points[0, 0, mi, oi]
+
+ # Loop over the time points, back calculating the R2eff values.
+ for j in range(num_points_i):
+
+ # If the number of spins are 1, do the fastest method by
dot product. Else do it by einstein summation.
+ if NS == 1:
+ # Set the spin index for one spin.
+ si = 0
# Extract the preformed matrix that rotate the
magnetization previous to spin-lock into the weff frame.
- M0_i= M0[0, si, mi, oi, j, :, 0]
-
- # This matrix is a propagator that will evolve the
magnetization with the matrix R.
- Rexpo_i = Rexpo_mat[0, si, mi, oi, j]
+ M0_i= M0_T[0, si, mi, oi, j]
# Extract from the pre-formed Magnetization evolution
matrix.
- Rexpo_M0_mat_i = Rexpo_M0_mat[0, si, mi, oi, j, :, 0]
+ Rexpo_M0_mat_i = Rexpo_M0_mat[0, si, mi, oi, j]
# Magnetization evolution.
- MA = dot(M0_i, Rexpo_M0_mat_i)
+ MA = dot(M0_i, Rexpo_M0_mat_i)[0, 0]
# The next lines calculate the R1rho using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
if MA <= 0.0 or isNaN(MA):
back_calc[0, si, mi, oi, j] = 1e99
else:
back_calc[0, si, mi, oi, j]= -inv_relax_time[0,
si, mi, oi, j] * log(MA)
+
+ # If there is multiple spin a clustered analysis.
+ else:
+ # Loop over spins.
+ for si in range(NS):
+ # Extract the preformed matrix that rotate the
magnetization previous to spin-lock into the weff frame.
+ M0_i= M0_T[0, si, mi, oi, j]
+
+ # Extract from the pre-formed Magnetization
evolution matrix.
+ Rexpo_M0_mat_i = Rexpo_M0_mat[0, si, mi, oi, j]
+
+ # Magnetization evolution.
+ MA = dot(M0_i, Rexpo_M0_mat_i)[0, 0]
+
+ MA_mat_i = MA_mat[0, si, mi, oi, j, 0, 0]
+
+ # Diff
+ diff = MA - MA_mat_i
+ if diff != 0.0:
+ print "oh no"
+
+ if MA_mat_i <= 0.0 or isNaN(MA_mat_i):
+ back_calc[0, si, mi, oi, j] = 1e99
+ else:
+ back_calc[0, si, mi, oi, j]=
-inv_relax_time[0, si, mi, oi, j] * log(MA_mat_i)
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-commits mailing list
relax-commits@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-commits
Re: r24292 - /branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py