r24233 - /branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_expanded.py -- June 21, 2014 - 10:21

Author: tlinnet
Date: Sat Jun 21 10:21:27 2014
New Revision: 24233

URL: http://svn.gna.org/viewcvs/relax?rev=24233&view=rev
Log:
Math domain fix for ns cpmg 2site expanded.

This is when t108 or t112 is zero, in the multidimensional array, a division 
error occurs.
The elements are first set to 1.0, to allow for computation.
Then elements are later replaced with 1e100.
Lastly, if the elements are not part of the "True" dispersion point 
structure, they are cleaned out.

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_cpmg_2site_expanded.py

Modified: branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_expanded.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_expanded.py?rev=24233&r1=24232&r2=24233&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_expanded.py   
(original)
+++ branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_expanded.py   Sat 
Jun 21 10:21:27 2014
@@ -236,7 +236,7 @@
 """
 
 # Python module imports.
-from numpy import exp, isfinite, fabs, power, log, min, sqrt, sum
+from numpy import any, exp, isfinite, fabs, power, log, min, sqrt, sum
 from numpy.ma import fix_invalid, masked_where
 
 
@@ -270,6 +270,8 @@
 
     # Flag to tell if values should be replaced if math function is violated.
     t_dw_zero = False
+    t_t108_zero = False
+    t_t112_zero = False
 
     # Catch parameter values that will result in no exchange, returning flat 
R2eff = R20 lines (when kex = 0.0, k_AB = 0.0).
     if pA == 1.0 or kex == 0.0:
@@ -348,7 +350,21 @@
     t99 = t92 + t96
     t102 = t99**2
     t108 = t62 * t88 + t82 * t31
+
+    # If t108 is zero.
+    mask_t108_zero = t108 == 0.0
+    if any(mask_t108_zero):
+        t_t108_zero = True
+        t108[mask_t108_zero] = 1.0
+
     t112 = sqrt(t98 - 2.0 * t99 * t97 + t102 + 2.0 * (t91 * t68 + t95 * t55) 
* t108)
+
+    # If t112 is zero.
+    mask_t112_zero = t112 == 0.0
+    if any(mask_t112_zero):
+        t_t112_zero = True
+        t112[mask_t112_zero] = 1.0
+
     t97_t99 = t97 + t99
     t97_nt99 = t97 - t99
     t113 = t97_nt99 - t112
@@ -376,6 +392,14 @@
     if t_dw_zero:
         back_calc[mask_dw_zero.mask] = r20[mask_dw_zero.mask]
 
+    # If t108 is zero.
+    if t_t108_zero:
+        back_calc[mask_t108_zero] = 1e100
+
+    # If t112 is zero.
+    if t_t112_zero:
+        back_calc[mask_t112_zero] = 1e100
+
     # 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)):