r23239 - in /branches/disp_speed: lib/dispersion/b14.py target_functions/relax_disp.py -- May 19, 2014 - 14:44

Author: tlinnet
Date: Mon May 19 14:44:10 2014
New Revision: 23239

URL: http://svn.gna.org/viewcvs/relax?rev=23239&view=rev
Log:
Math-domain catching for model B14.

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

This is to implement catching of math domain errors, before they occur.
These can be found via the --numpy-raise function to the systemtests.

To make the code look clean, the class object "back_calc" is no longer
being updated per time point, but is updated in the relax_disp target 
function in
one go.

Modified:
    branches/disp_speed/lib/dispersion/b14.py
    branches/disp_speed/target_functions/relax_disp.py

Modified: branches/disp_speed/lib/dispersion/b14.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/b14.py?rev=23239&r1=23238&r2=23239&view=diff
==============================================================================
--- branches/disp_speed/lib/dispersion/b14.py   (original)
+++ branches/disp_speed/lib/dispersion/b14.py   Mon May 19 14:44:10 2014
@@ -110,12 +110,12 @@
 """
 
 # Python module imports.
-from numpy import arccosh, array, cos, cosh, isfinite, log, power, sin, 
sinh, sqrt, sum
+from numpy import abs, arccosh, array, cos, cosh, isfinite, log, min, power, 
sin, sinh, sqrt, sum
 
 # Repetitive calculations (to speed up calculations).
 g_fact = 1/sqrt(2)
 
-def r2eff_B14(r20a=None, r20b=None, pA=None, pB=None, dw=None, kex=None, 
k_AB=None, k_BA=None, ncyc=None, inv_tcpmg=None, tcp=None, back_calc=None, 
num_points=None):
+def r2eff_B14(r20a=None, r20b=None, pA=None, pB=None, dw=None, kex=None, 
k_AB=None, k_BA=None, ncyc=None, inv_tcpmg=None, tcp=None, num_points=None):
     """Calculate the R2eff values for the CR72 model.
 
     See the module docstring for details.
@@ -143,9 +143,7 @@
     @type inv_tcpmg:        float
     @keyword tcp:           The tau_CPMG times (1 / 4.nu1).
     @type tcp:              numpy rank-1 float array
-    @keyword back_calc:     The array for holding the back calculated R2eff 
values.  Each element corresponds to one of the CPMG nu1 frequencies.
-    @type back_calc:        numpy rank-1 float array
-    @keyword num_points:    The number of points on the dispersion curve, 
equal to the length of the cpmg_frqs and back_calc arguments.
+    @keyword num_points:    The number of points on the dispersion curve, 
equal to the length of the cpmg_frqs.
     @type num_points:       int
     """
 
@@ -215,12 +213,24 @@
     # Real. The v_1c in paper.
     v1c = F0 * cosh(E0) - F2 * cos(E2)
 
+    # Catch math domain error of sqrt(neg val).
+    # This is when abs(v1c) =< 1.
+    if min(abs(v1c)) <= 1.:
+        R2eff = array([1e100]*num_points)
+        return 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 )
+
+    # Catch math domain error of log(neg val or zero).
+    # This is when Tog.real =< 0.
+    if min(Tog.real) <= 0.:
+        R2eff = array([1e100]*num_points)
+        return R2eff
 
     ## -1/Trel * log(LpreDyn).
     # Rpre = (r20a + r20b + kex) / 2.0
@@ -240,6 +250,4 @@
     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]
+    return R2eff

Modified: branches/disp_speed/target_functions/relax_disp.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_speed/target_functions/relax_disp.py?rev=23239&r1=23238&r2=23239&view=diff
==============================================================================
--- branches/disp_speed/target_functions/relax_disp.py  (original)
+++ branches/disp_speed/target_functions/relax_disp.py  Mon May 19 14:44:10 
2014
@@ -439,7 +439,7 @@
                         aliased_dw = dw_frq
 
                     # Back calculate the R2eff values.
-                    r2eff_B14(r20a=R20A[r20_index], r20b=R20B[r20_index], 
pA=pA, pB=pB, dw=aliased_dw, kex=kex, k_AB=k_AB, k_BA=k_BA, 
ncyc=self.power[ei][mi], inv_tcpmg=self.inv_relax_times[ei][mi], 
tcp=self.tau_cpmg[ei][mi], back_calc=self.back_calc[ei][si][mi][0], 
num_points=self.num_disp_points[ei][si][mi][0])
+                    self.back_calc[ei][si][mi][0] = 
r2eff_B14(r20a=R20A[r20_index], r20b=R20B[r20_index], pA=pA, pB=pB, 
dw=aliased_dw, kex=kex, k_AB=k_AB, k_BA=k_BA, ncyc=self.power[ei][mi], 
inv_tcpmg=self.inv_relax_times[ei][mi], tcp=self.tau_cpmg[ei][mi], 
num_points=self.num_disp_points[ei][si][mi][0])
 
                     # For all missing data points, set the back-calculated 
value to the measured values so that it has no effect on the chi-squared 
value.
                     for di in range(self.num_disp_points[ei][si][mi][0]):