Author: tlinnet
Date: Mon May 19 03:20:46 2014
New Revision: 23223
URL: http://svn.gna.org/viewcvs/relax?rev=23223&view=rev
Log:
Speed-up of model MP05.
task #7793: (https://gna.org/task/?7793) Speed-up of dispersion models.
The change in systemtest is:
test_tp02_data_to_mp05
10.750s -> 6.644s
Modified:
branches/disp_speed/lib/dispersion/mp05.py
Modified: branches/disp_speed/lib/dispersion/mp05.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/mp05.py?rev=23223&r1=23222&r2=23223&view=diff
==============================================================================
--- branches/disp_speed/lib/dispersion/mp05.py (original)
+++ branches/disp_speed/lib/dispersion/mp05.py Mon May 19 03:20:46 2014
@@ -60,7 +60,7 @@
"""
# Python module imports.
-from math import atan2, sin
+from numpy import arctan2, isfinite, sin, sum
def r1rho_MP05(r1rho_prime=None, omega=None, offset=None, pA=None, pB=None,
dw=None, kex=None, R1=0.0, spin_lock_fields=None, spin_lock_fields2=None,
back_calc=None, num_points=None):
@@ -104,39 +104,37 @@
phi_ex = pA * pB * dw**2
numer = phi_ex * kex
- # Loop over the dispersion points, back calculating the R1rho values.
+ # We assume that A resonates at 0 [s^-1], without loss of generality.
+ W = pA*Wa + pB*Wb # Pop-averaged Larmor
frequency [s^-1].
+ da = Wa - offset # Offset of spin-lock from A.
+ db = Wb - offset # Offset of spin-lock from B.
+ d = W - offset # Offset of spin-lock from
pop-average.
+ waeff2 = spin_lock_fields2 + da**2 # Effective field at A.
+ wbeff2 = spin_lock_fields2 + db**2 # Effective field at B.
+ weff2 = spin_lock_fields2 + d**2 # Effective field at
pop-average.
+
+ # The rotating frame flip angle.
+ theta = arctan2(spin_lock_fields, d)
+
+ # Repetitive calculations (to speed up calculations).
+ sin_theta2 = sin(theta)**2
+ R1_cos_theta2 = R1 * (1.0 - sin_theta2)
+ R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2
+
+ # Denominator.
+ waeff2_wbeff2 = waeff2*wbeff2
+ fact = 1.0 + 2.0*kex2*(pA*waeff2 + pB*wbeff2) / (waeff2_wbeff2 +
weff2*kex2)
+ denom = waeff2_wbeff2/weff2 + kex2 - sin_theta2*phi_ex*(fact)
+
+ # R1rho calculation.
+ R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer /
denom
+
+ # 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(R1rho)):
+ R1rho = array([1e100]*num_points)
+
+ # Parse back the value to update the back_calc class object.
for i in range(num_points):
- # We assume that A resonates at 0 [s^-1], without loss of generality.
- W = pA*Wa + pB*Wb # Pop-averaged Larmor
frequency [s^-1].
- da = Wa - offset # Offset of spin-lock
from A.
- db = Wb - offset # Offset of spin-lock
from B.
- d = W - offset # Offset of spin-lock
from pop-average.
- waeff2 = spin_lock_fields2[i] + da**2 # Effective field at A.
- wbeff2 = spin_lock_fields2[i] + db**2 # Effective field at B.
- weff2 = spin_lock_fields2[i] + d**2 # Effective field at
pop-average.
+ back_calc[i] = R1rho[i]
- # The rotating frame flip angle.
- theta = atan2(spin_lock_fields[i], d)
-
- # Repetitive calculations (to speed up calculations).
- sin_theta2 = sin(theta)**2
- R1_cos_theta2 = R1 * (1.0 - sin_theta2)
- R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2
-
- # Catch zeros (to avoid pointless mathematical operations).
- if numer == 0.0:
- back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2
- continue
-
- # Denominator.
- waeff2_wbeff2 = waeff2*wbeff2
- fact = 1.0 + 2.0*kex2*(pA*waeff2 + pB*wbeff2) / (waeff2_wbeff2 +
weff2*kex2)
- denom = waeff2_wbeff2/weff2 + kex2 - sin_theta2*phi_ex*(fact)
-
- # Avoid divide by zero.
- if denom == 0.0:
- back_calc[i] = 1e100
- continue
-
- # R1rho calculation.
- back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 *
numer / denom
r23223 - /branches/disp_speed/lib/dispersion/mp05.py