r23220 - /branches/disp_speed/lib/dispersion/b14.py -- May 19, 2014 - 00:51

Author: tlinnet
Date: Mon May 19 00:51:20 2014
New Revision: 23220

URL: http://svn.gna.org/viewcvs/relax?rev=23220&view=rev
Log:
Huge speed-up of model B14.

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

Time test for systemtests:

test_baldwin_synthetic
2.626s -> 1.990s

test_baldwin_synthetic_full
18.326s -> 13.742s

This is won by not checking single values in the R2eff 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/b14.py

Modified: branches/disp_speed/lib/dispersion/b14.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/b14.py?rev=23220&r1=23219&r2=23220&view=diff
==============================================================================
--- branches/disp_speed/lib/dispersion/b14.py   (original)
+++ branches/disp_speed/lib/dispersion/b14.py   Mon May 19 00:51:20 2014
@@ -110,8 +110,7 @@
 """
 
 # Python module imports.
-import numpy
-from numpy import arccosh, array, cos, cosh, in1d, log, nonzero, sin, sinh, 
sqrt, power
+from numpy import arccosh, array, cos, cosh, isfinite, log, power, sin, 
sinh, sqrt, sum
 
 # Repetitive calculations (to speed up calculations).
 g_fact = 1/sqrt(2)
@@ -216,70 +215,31 @@
     # Real. The v_1c in paper.
     v1c = F0 * cosh(E0) - F2 * cos(E2)
 
-    # Catch devision with zero in y, when v3 = 0. v3 is 0, when v1c = 1.
-    # If no 1.0, perform normally.
-    if not in1d(1.0, v1c):
-        # Exact result for v2v3.
-        v3 = sqrt(v1c**2 - 1.)
-
-        y = power( (v1c - v3) / (v1c + v3), ncyc)
-
-        Tog = 0.5 * (1. + y) + (1. - y) * v5 / (2. * v3 * N )
-
-        # Find where Tog has negative values.
-        neg_index = nonzero(Tog.real < 0.0)[0]
-
-        # Do normal calculation
-        if len(neg_index) == 0:
-            ## -1/Trel * log(LpreDyn).
-            # Rpre = (r20a + r20b + kex) / 2.0
-
-            ## Carver and Richards (1972)
-            # R2eff_CR72 = Rpre - inv_tcpmg * ncyc *  arccosh(v1c.real)
-
-            ## Baldwin final.
-            # Estimate R2eff. relax_time = Trel = 1/inv_tcpmg.
-            # R2eff = R2eff_CR72 - inv_tcpmg * log(Tog.real)
-
-            # Fastest calculation.
-            R2eff = (r20a + r20b + kex) / 2.0  - inv_tcpmg * ( ncyc *  
arccosh(v1c.real) + log(Tog.real) )
-
-            # Loop over the time points, back calculating the R2eff values.
-            for i in range(num_points):
-
-                # Put values back.
-                back_calc[i] = R2eff[i]
-
-        else:
-            # Loop over each point.
-            for i in range(num_points):
-
-                # Return large value
-                if i in neg_index:
-                    back_calc[i] = 1e100        
-
-                else:
-                    v3 = sqrt(v1c[i]**2 - 1.)
-                    y = power( (v1c[i] - v3) / (v1c[i] + v3), ncyc[i])
-                    Tog = 0.5 * (1. + y) + (1. - y) * v5[i] / (2. * v3 * N )
-                    R2eff = (r20a + r20b + kex) / 2.0  - inv_tcpmg * ( 
ncyc[i] *  arccosh(v1c[i].real) + log(Tog.real) )
-                    back_calc[i] = R2eff
-
-    # This section is for catching math domain errors.
-    else:
-        # Find index where 
-        one_indexes = nonzero(v1c == 1.0)[0]
-
-        # Loop over each point.
-        for i in range(num_points):
-
-            # Return large value
-            if i in one_indexes:
-                back_calc[i] = 1e100
-
-            else:
-                v3 = sqrt(v1c[i]**2 - 1.)
-                y = power( (v1c[i] - v3) / (v1c[i] + v3), ncyc[i])
-                Tog = 0.5 * (1. + y) + (1. - y) * v5[i] / (2. * v3 * N )
-                R2eff = (r20a + r20b + kex) / 2.0  - inv_tcpmg * ( ncyc[i] * 
 arccosh(v1c[i].real) + log(Tog.real) )
-                back_calc[i] = R2eff
+    # Exact result for v2v3.
+    v3 = sqrt(v1c**2 - 1.)
+
+    y = power( (v1c - v3) / (v1c + v3), ncyc)
+
+    Tog = 0.5 * (1. + y) + (1. - y) * v5 / (2. * v3 * N )
+
+    ## -1/Trel * log(LpreDyn).
+    # Rpre = (r20a + r20b + kex) / 2.0
+
+    ## Carver and Richards (1972)
+    # R2eff_CR72 = Rpre - inv_tcpmg * ncyc *  arccosh(v1c.real)
+
+    ## Baldwin final.
+    # Estimate R2eff. relax_time = Trel = 1/inv_tcpmg.
+    # R2eff = R2eff_CR72 - inv_tcpmg * log(Tog.real)
+
+    # Fastest calculation.
+    R2eff = (r20a + r20b + kex) / 2.0  - inv_tcpmg * ( ncyc *  
arccosh(v1c.real) + log(Tog.real) )
+
+    # 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(R2eff)):
+        R2eff = array([1e100]*num_points)
+
+    # Parse back the value to update the back_calc class object.
+    for i in range(num_points):
+        back_calc[i] = R2eff[i]