Re: r23221 - /branches/disp_speed/lib/dispersion/tp02.py -- May 19, 2014

Hi Troels,

Here again you can avoid the divide by zero problems by having:

"""
from numpy import max, abs

if min(abs(denom)) == 0:
    back_calc[:] = array([1e100]*num_points)
"""

Then these tests should then pass when running with --numpy-raise.
You should also not see any significant speed changes (at least I
didn't with a quick test).

Regards,

Edward



On 19 May 2014 01:19,  <tlinnet@xxxxxxxxxxxxx> wrote:

Author: tlinnet
Date: Mon May 19 01:19:02 2014
New Revision: 23221

URL: http://svn.gna.org/viewcvs/relax?rev=23221&view=rev
Log:
Speed-up of model TP02.

task #7793: (https://gna.org/task/?7793) Speed-up of dispersion models.

The change for running systemtest is:
test_curve_type_r1rho_fixed_time
0.057s -> 0.049s

test_tp02_data_to_ns_r1rho_2site
10.539s -> 10.456s

test_tp02_data_to_tp02
8.608s -> 5.727s

This is won by not checking single values in the R1rho array for math domain
errors, but calculating all steps, and in one single round check for finite 
values.
If just one non-finite value is found, the whole array is returned with a 
large
penalty of 1e100.

This makes all calculations be the fastest numpy array way.

Modified:
    branches/disp_speed/lib/dispersion/tp02.py

Modified: branches/disp_speed/lib/dispersion/tp02.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/tp02.py?rev=23221&r1=23220&r2=23221&view=diff
==============================================================================
--- branches/disp_speed/lib/dispersion/tp02.py  (original)
+++ branches/disp_speed/lib/dispersion/tp02.py  Mon May 19 01:19:02 2014
@@ -60,7 +60,7 @@
 """

 # Python module imports.
-from math import atan2, sin
+from numpy import arctan2, isfinite, sin, sum


 def r1rho_TP02(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):
@@ -110,34 +110,31 @@
     # The numerator.
     numer = pA * pB * dw**2 * kex

-    # Loop over the dispersion points, back calculating the R1rho values.
+    # We assume that A resonates at 0 [s^-1], without loss of generality.
+    waeff2 = spin_lock_fields2 + da2       # Effective field at A.
+    wbeff2 = spin_lock_fields2 + db2       # Effective field at B.
+    weff2 = spin_lock_fields2 + d2         # 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.
+    denom = waeff2 * wbeff2 / weff2 + kex2
+    #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2
+
+    # 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.
-        waeff2 = spin_lock_fields2[i] + da2       # Effective field at A.
-        wbeff2 = spin_lock_fields2[i] + db2       # Effective field at B.
-        weff2 = spin_lock_fields2[i] + d2         # Effective field at 
pop-average.
-
-        # 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.
-        denom = waeff2 * wbeff2 / weff2 + kex2
-        #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2
-
-        # 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
+        back_calc[i] = R1rho[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: r23221 - /branches/disp_speed/lib/dispersion/tp02.py

Header

Content

Related Messages