Author: tlinnet
Date: Wed Jun 25 02:14:36 2014
New Revision: 24294
URL: http://svn.gna.org/viewcvs/relax?rev=24294&view=rev
Log:
Speeded up the code of NS r1rho 2site.
This was essential done to numpy einsum, and doing the dot operations in
multiple dimensions.
It was though necessary to realize, that to do the proper dot product
operations, the outer two
axis if M0 should be swapped, by rolling the outer axis one back.
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
branches/disp_spin_speed/target_functions/relax_disp.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=24294&r1=24293&r2=24294&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:36 2014
@@ -50,10 +50,10 @@
"""
# Python module imports.
-from numpy import array, cos, dot, einsum, float64, log, multiply, rollaxis,
sin, sum
+from numpy import array, einsum, float64, isfinite, log, min, multiply,
rollaxis, sin, sum
+from numpy.ma import fix_invalid, masked_less
# relax module imports.
-from lib.float import isNaN
from lib.dispersion.matrix_exponential import
matrix_exponential_rank_NE_NS_NM_NO_ND_x_x
# Repetitive calculations (to speed up calculations).
@@ -187,14 +187,16 @@
return matrix
-def ns_r1rho_2site(M0=None, r1rho_prime=None, omega=None, offset=None,
r1=0.0, pA=None, dw=None, kex=None, spin_lock_fields=None, relax_time=None,
inv_relax_time=None, back_calc=None, num_points=None):
+def ns_r1rho_2site(M0=None, M0_T=None, r1rho_prime=None, omega=None,
offset=None, r1=0.0, pA=None, dw=None, kex=None, spin_lock_fields=None,
relax_time=None, inv_relax_time=None, back_calc=None, num_points=None):
"""The 2-site numerical solution to the Bloch-McConnell equation for
R1rho data.
This function calculates and stores the R1rho values.
@keyword M0: This is a vector that contains the initial
magnetizations corresponding to the A and B state transverse magnetizations.
- @type M0: numpy float64, rank-1, 7D array
+ @type M0: numpy float array of rank
[NE][NS][NM][NO][ND][6][1]
+ @keyword M0_T: This is a vector that contains the initial
magnetizations corresponding to the A and B state transverse magnetizations,
where the outer two axis has been swapped for efficient dot operations.
+ @type M0_T: numpy float array of rank
[NE][NS][NM][NO][ND][1][6]
@keyword r1rho_prime: The R1rho_prime parameter value (R1rho with
no exchange).
@type r1rho_prime: numpy float array of rank
[NE][NS][NM][NO][ND]
@keyword omega: The chemical shift for the spin in rad/s.
@@ -238,43 +240,22 @@
# Magnetization evolution.
Rexpo_M0_mat = einsum('...ij,...jk', Rexpo_mat, M0)
- # Transpose M0, to prepare for dot operation. Roll the last axis one
back.
- M0_T = rollaxis(M0, 6, 5)
-
- if NS == 1:
- # 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.
- # 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_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]
-
- # 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)
-
- else:
- # 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)
-
-
+ # Magnetization evolution, which include all dimensions.
+ MA_mat = einsum('...ij,...jk', M0_T, Rexpo_M0_mat)[:, :, :, :, :, 0, 0]
+
+ # Insert safe checks.
+ if min(MA_mat) < 0.0:
+ mask_min_MA_mat = masked_less(MA_mat, 0.0)
+ # Fill with high values.
+ MA_mat[mask_min_MA_mat.mask] = 1e100
+
+ # Do back calculation.
+ back_calc[:] = -inv_relax_time * log(MA_mat)
+
+ # Catch errors, taking a sum over array is the fastest way to check for
+ # +/- inf (infinity) and nan (not a number).
+ if not isfinite(sum(back_calc)):
+ # Replaces nan, inf, etc. with fill value.
+ fix_invalid(back_calc, copy=False, fill_value=1e100)
+
+
Modified: branches/disp_spin_speed/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/target_functions/relax_disp.py?rev=24294&r1=24293&r2=24294&view=diff
==============================================================================
--- branches/disp_spin_speed/target_functions/relax_disp.py (original)
+++ branches/disp_spin_speed/target_functions/relax_disp.py Wed Jun 25
02:14:36 2014
@@ -26,8 +26,7 @@
# Python module imports.
from copy import deepcopy
-from math import pi
-from numpy import add, array, arctan2, asarray, cos, complex64, dot,
float64, int16, max, multiply, ones, sin, sqrt, sum, tile, zeros
+from numpy import arctan2, cos, dot, float64, int16, multiply, ones,
rollaxis, pi, sin, sum, zeros
from numpy.ma import masked_equal
# relax module imports.
@@ -415,6 +414,9 @@
M0_2[2, 0] = 1
M0_cos = multiply.outer( cos(theta_mat), M0_2 )
self.M0 = M0_sin + M0_cos
+ # Transpose M0, to prepare for dot operation. Roll the last axis
one back, corresponds to a transpose for the outer two axis.
+ self.M0_T = rollaxis(self.M0, 6, 5)
+
if model in [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]:
self.M0 = zeros(9, float64)
# Offset of spin-lock from A.
@@ -1598,7 +1600,7 @@
self.r20_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE,
self.NS, self.NM), self.no_nd_ones )
# Back calculate the R2eff values.
- ns_r1rho_2site(M0=self.M0, r1rho_prime=self.r20_struct,
omega=self.chemical_shifts, offset=self.offset, r1=self.r1, pA=pA,
dw=self.dw_struct, kex=kex, spin_lock_fields=self.spin_lock_omega1,
relax_time=self.relax_times, inv_relax_time=self.inv_relax_times,
back_calc=self.back_calc, num_points=self.num_disp_points)
+ ns_r1rho_2site(M0=self.M0, M0_T=self.M0_T,
r1rho_prime=self.r20_struct, omega=self.chemical_shifts, offset=self.offset,
r1=self.r1, pA=pA, dw=self.dw_struct, kex=kex,
spin_lock_fields=self.spin_lock_omega1, relax_time=self.relax_times,
inv_relax_time=self.inv_relax_times, back_calc=self.back_calc,
num_points=self.num_disp_points)
# Clean the data for all values, which is left over at the end of
arrays.
self.back_calc = self.back_calc*self.disp_struct
r24294 - in /branches/disp_spin_speed: lib/dispersion/ns_r1rho_2site.py target_functions/relax_disp.py