Re: r23091 - /trunk/lib/dispersion/b14.py -- May 08, 2014

Very much so :)  On a related note, you would get insane speed ups of
the dispersion model if one or more of the loops in the target
functions could be brought into these lib.dispersion modules and then
this exact logic used.  However I have never managed to make this work
- though I haven't tried too hard.  It would require quite large numpy
data structures with multiply replicated data.  Such a revolutionary
idea would really need its own real subversion branch, and a lot of
playing around testing ideas.

Regards,

Edward



On 8 May 2014 19:03, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:

Yeah.

I tried really hard to stay in numpy land as long as possible.

Are you not here introducing the rather stupid mistake of always
calculating in a loop?

I have seen the same for CR72.

Best
Troels

2014-05-08 18:50 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:

Hi,

For a test, and as a reference, I replaced all lines after the v5 =  line 
with:

"""
    # Real. The v_1c in paper.
    v1c = F0 * cosh(E0) - F2 * cos(E2)

    # Loop over each point.
    for i in range(num_points):
        # Catch devision with zero in y, when v3 = 0. v3 is 0, when v1c = 
1.
        # If no 1.0, perform normally.
        if v1c[i] == 1:
            back_calc[i] = 1e100
            continue

        # Exact result for v2v3.
        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 )

        # Find where Tog has negative values.
        if Tog.real < 0.0:
            back_calc[i] = 1e100
            continue

        # Fastest calculation.
        back_calc[i] = (r20a + r20b + kex) / 2.0  - inv_tcpmg * (
ncyc[i] *  arccosh(v1c[i].real) + log(Tog.real) )
"""

This is far more compact and easier to read.  However the
Relax_disp.test_baldwin_synthetic_full system test time on my system
went from 7.648s to 11.155s!  This difference is because the old code
uses the numpy speed advantage of operating on numpy arrays, and the
above code does not.

Regards,

Edward

On 8 May 2014 18:35,  <tlinnet@xxxxxxxxxxxxx> wrote:

Author: tlinnet
Date: Thu May  8 18:35:51 2014
New Revision: 23091

URL: http://svn.gna.org/viewcvs/relax?rev=23091&view=rev
Log:
Made solutions for math domain error.

Prevented to take log of negative values, and division by zero.

This though slows the implementation down.

Systemtest Relax_disp.test_baldwin_synthetic_full went from 6.x seconds 
to 8-9.x seconds.

sr #3154: (https://gna.org/support/?3154) Implementation of Baldwin 
(2014) B14 model - 2-site exact solution model for all time scales.

This follows the tutorial for adding relaxation dispersion models at:
http://wiki.nmr-relax.com/Tutorial_for_adding_relaxation_dispersion_models_to_relax#Debugging

Modified:
    trunk/lib/dispersion/b14.py

Modified: trunk/lib/dispersion/b14.py
URL: 
http://svn.gna.org/viewcvs/relax/trunk/lib/dispersion/b14.py?rev=23091&r1=23090&r2=23091&view=diff
==============================================================================
--- trunk/lib/dispersion/b14.py (original)
+++ trunk/lib/dispersion/b14.py Thu May  8 18:35:51 2014
@@ -111,7 +111,7 @@

 # Python module imports.
 import numpy
-from numpy import arccosh, cos, cosh, log, sin, sinh, sqrt, power
+from numpy import arccosh, array, cos, cosh, in1d, log, nonzero, sin, 
sinh, sqrt, power

 # Repetitive calculations (to speed up calculations).
 g_fact = 1/sqrt(2)
@@ -201,43 +201,85 @@
     # Mixed term (complex) (E0 - iE2)/2.
     E1 = (g3 - g4*1j) * tcp

+    # Complex.
+    v1s = F0 * sinh(E0) - F2 * sin(E2)*1j
+
+    # -2 * oG * t2.
+    v4 = F1b * (-alpha_m - g3 ) + F1b * (dw - g4)*1j
+
+    # Complex.
+    ex1c = sinh(E1)
+
+    # Off diagonal common factor. sinh fuctions.
+    v5 = (-deltaR2 + kex + dw*1j) * v1s - 2. * (v4 + k_AB * F1a_plus_b) 
* ex1c
+
     # Real. The v_1c in paper.
     v1c = F0 * cosh(E0) - F2 * cos(E2)

-    # Complex.
-    v1s = F0 * sinh(E0) - F2 * sin(E2)*1j
-
-    # Exact result for v2v3.
-    v3 = sqrt(v1c**2 - 1.)
-
-    # -2 * oG * t2.
-    v4 = F1b * (-alpha_m - g3 ) + F1b * (dw - g4)*1j
-
-    # Complex.
-    ex1c = sinh(E1)
-
-    # Off diagonal common factor. sinh fuctions.
-    v5 = (-deltaR2 + kex + dw*1j) * v1s - 2. * (v4 + k_AB * F1a_plus_b) 
* ex1c
-
-    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) )
-
-    # Loop over the time points, back calculating the R2eff values.
-    for i in range(num_points):
-
-        # The full formula.
-        back_calc[i] = R2eff[i]
+    # 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


_______________________________________________
relax (http://www.nmr-relax.com)

This is the relax-commits mailing list
relax-commits@xxxxxxx

To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-commits


_______________________________________________
relax (http://www.nmr-relax.com)

This is the relax-devel mailing list
relax-devel@xxxxxxx

To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-devel

Re: r23091 - /trunk/lib/dispersion/b14.py

Header

Content

Related Messages