r23937 - /branches/disp_spin_speed/lib/dispersion/it99.py -- June 13, 2014 - 17:31

Author: tlinnet
Date: Fri Jun 13 17:31:31 2014
New Revision: 23937

URL: http://svn.gna.org/viewcvs/relax?rev=23937&view=rev
Log:
Methods to replace math domain errors in model IT99, has been replaced with 
numpy masks.

number of points has been removed, as the masks utility replaces this.
pB is now moved to be calculated inside. This makes the lib function simpler.

Documentation is also fixed.

Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion 
models for Clustered analysis.

Modified:
    branches/disp_spin_speed/lib/dispersion/it99.py

Modified: branches/disp_spin_speed/lib/dispersion/it99.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/it99.py?rev=23937&r1=23936&r2=23937&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/it99.py     (original)
+++ branches/disp_spin_speed/lib/dispersion/it99.py     Fri Jun 13 17:31:31 
2014
@@ -73,33 +73,47 @@
 """
 
 # Python module imports.
-from numpy import array, isfinite, sqrt, sum
+from numpy import array, isfinite, fabs, sqrt, sum
+from numpy.ma import fix_invalid, masked_where
 
 
-def r2eff_IT99(r20=None, pA=None, pB=None, dw=None, tex=None, 
cpmg_frqs=None, back_calc=None, num_points=None):
+def r2eff_IT99(r20=None, pA=None, dw=None, dw_orig=None, tex=None, 
cpmg_frqs=None, back_calc=None):
     """Calculate the R2eff values for the IT99 model.
 
     See the module docstring for details.
 
 
     @keyword r20:           The R20 parameter value (R2 with no exchange).
-    @type r20:              float
+    @type r20:              numpy float array of rank [NE][NS][[NM][NO][ND]
     @keyword pA:            The population of state A.
     @type pA:               float
-    @keyword pB:            The population of state B.
-    @type pB:               float
     @keyword dw:            The chemical exchange difference between states 
A and B in rad/s.
-    @type dw:               float
+    @type dw:               numpy float array of rank [NE][NS][[NM][NO][ND]
+    @keyword dw_orig:       The chemical exchange difference between states 
A and B in ppm. This is only for faster checking of zero value, which result 
in no exchange.
+    @type dw_orig:          numpy float array of rank-1
     @keyword tex:           The tex parameter value (the time of exchange in 
s/rad).
     @type tex:              float
     @keyword cpmg_frqs:     The CPMG nu1 frequencies.
-    @type cpmg_frqs:        numpy rank-1 float array
-    @type tcp:              numpy rank-1 float array
+    @type cpmg_frqs:        numpy float array of rank [NE][NS][[NM][NO][ND]
     @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.
-    @type num_points:       int
+    @type back_calc:        numpy float array of rank [NE][NS][[NM][NO][ND]
     """
+
+    # Flag to tell if values should be replaced if numer is zero.
+    t_dw_zero = False
+
+    # Catch divide with zeros (to avoid pointless mathematical operations).
+    if tex == 0.0 or pA == 1.0:
+        back_calc[:] = r20
+        return
+
+    # Test if dw is zero. Wait for replacement, since this is spin specific.
+    if min(fabs(dw_orig)) == 0.0:
+        t_dw_zero = True
+        mask_dw_zero = masked_where(dw == 0.0, dw)
+
+    # Parameter conversions.
+    pB = 1.0 - pA
 
     # Repetitive calculations (to speed up calculations).
     dw2 = dw**2
@@ -110,12 +124,6 @@
     # The numerator.
     numer = padw2 * pB * tex
 
-    # Catch zeros (to avoid pointless mathematical operations).
-    # This will result in no exchange, returning flat lines.
-    if numer == 0.0:
-        back_calc[:] = array([r20]*num_points)
-        return
-
     # The effective rotating frame field.
     omega_1eff4 = 2304.0 * cpmg_frqs**4
 
@@ -124,11 +132,15 @@
     denom = 1.0 + omega_a2 * tex2
 
     # R2eff calculation.
-    R2eff = r20 + numer / denom
+    back_calc[:] = r20 + numer / denom
+
+    # Replace data in array.
+    # If dw is zero.
+    if t_dw_zero:
+        back_calc[mask_dw_zero.mask] = r20[mask_dw_zero.mask]
 
     # 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)
-
-    back_calc[:] = R2eff
+    if not isfinite(sum(back_calc)):
+        # Replaces nan, inf, etc. with fill value.
+        fix_invalid(back_calc, copy=False, fill_value=1e100)