r23797 - /branches/disp_spin_speed/lib/dispersion/cr72.py -- June 10, 2014 - 18:41

Author: tlinnet
Date: Tue Jun 10 18:41:08 2014
New Revision: 23797

URL: http://svn.gna.org/viewcvs/relax?rev=23797&view=rev
Log:
Important fix for replacing values if eta_pos > 700 is violated.

This fixes systemtest: Relax_disp.test_sod1wt_t25_to_cr72, which failed after 
making kex to a numpy float.

The trick is to make a numpy mask which stores the position where to replace 
the values.
Then replace the values just before last return.

This makes sure, that not all values are changed.

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

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

Modified: branches/disp_spin_speed/lib/dispersion/cr72.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/cr72.py?rev=23797&r1=23796&r2=23797&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/cr72.py     (original)
+++ branches/disp_spin_speed/lib/dispersion/cr72.py     Tue Jun 10 18:41:08 
2014
@@ -93,7 +93,7 @@
 
 # Python module imports.
 from numpy import allclose, arccosh, array, cos, cosh, isfinite, isnan, min, 
max, ndarray, ones, sqrt, sum, zeros
-import numpy.ma as ma
+from numpy.ma import masked_greater_equal
 
 # Repetitive calculations (to speed up calculations).
 eta_scale = 2.0**(-3.0/2.0)
@@ -127,6 +127,9 @@
     #mask_pA_t = False
     #mask_kex_t = False
 
+    # Flag to tell if values should be replaced if max_etapos in cosh 
function is violated.
+    t_max_etapos = False
+
     # Determine if calculating in numpy rank-1 float array, of higher 
dimensions.
     rank_1 = True
     if isinstance(num_points, ndarray):
@@ -193,10 +196,15 @@
         if rank_1:
             back_calc[:] = array([r20a]*num_points)
             return
-        # For higher dimensions, return same structure.
+        # For higher dimensions, find the mask to replace values.
+        # Reset to 1.0 and wait for replacement to later.
         else:
-            back_calc[:] = r20a
-            return
+            # Set the flag to tell to replace values.
+            t_max_etapos = True
+            # Find the mask, where to replace values.            
+            mask_max_etapos = masked_greater_equal(etapos, 700.0)
+            # To prevent math errors, set etapos to 1.
+            etapos[mask_max_etapos.mask] = 1.0
 
     # The arccosh argument - catch invalid values.
     fact = Dpos * cosh(etapos) - Dneg * cos(etaneg)
@@ -218,6 +226,8 @@
     #    R2eff[mask_pA] = r20a
     #if mask_kex_t:
     #    R2eff[mask_kex] = r20a
+    if t_max_etapos:
+        R2eff[mask_max_etapos.mask] = r20a[mask_max_etapos.mask]
 
     # Catch errors, taking a sum over array is the fastest way to check for
     # +/- inf (infinity) and nan (not a number).