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

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