Re: r23224 - /branches/disp_speed/lib/dispersion/mmq_cr72.py -- May 19, 2014 - 10:42

Hi Troels,

This one will be much harder to properly catch all numpy errors.
However you could try something like:

"""
from numpy import abs, max

if max(abs(etapos)) > 200:
    back_calc[:] = array([1e100]*num_points)
    return
"""

This might be sufficient to replace your isfinite() block.

Regards,

Edward



On 19 May 2014 03:20,  <tlinnet@xxxxxxxxxxxxx> wrote:
> Author: tlinnet
> Date: Mon May 19 03:20:47 2014
> New Revision: 23224
>
> URL: http://svn.gna.org/viewcvs/relax?rev=23224&view=rev
> Log:
> Speed-up of model MMQ CR72.
>
> task #7793: (https://gna.org/task/?7793) Speed-up of dispersion models.
>
> Change in systemtest:
> test_sprangers_data_to_mmq_cr72
> 9.892s ->  4.121s
>
> Modified:
>     branches/disp_speed/lib/dispersion/mmq_cr72.py
>
> Modified: branches/disp_speed/lib/dispersion/mmq_cr72.py
> URL: 
> http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/mmq_cr72.py?rev=23224&r1=23223&r2=23224&view=diff
> ==============================================================================
> --- branches/disp_speed/lib/dispersion/mmq_cr72.py      (original)
> +++ branches/disp_speed/lib/dispersion/mmq_cr72.py      Mon May 19 03:20:47 
> 2014
> @@ -47,7 +47,7 @@
>  """
>
>  # Python module imports.
> -from numpy import arccosh, cos, cosh, log, sin, sqrt
> +from numpy import arccosh, cos, cosh, isfinite, log, sin, sqrt, sum
>
>
>  def r2eff_mmq_cr72(r20=None, pA=None, pB=None, dw=None, dwH=None, 
> kex=None, k_AB=None, k_BA=None, cpmg_frqs=None, inv_tcpmg=None, tcp=None, 
> back_calc=None, num_points=None, power=None):
> @@ -122,27 +122,31 @@
>      etapos_part = eta_scale * sqrt(Psi + sqrt_psi2_zeta2)
>      etaneg_part = eta_scale * sqrt(-Psi + sqrt_psi2_zeta2)
>
> -    # Loop over the time points, back calculating the R2eff values.
> +    # The full eta+/- values.
> +    etapos = etapos_part / cpmg_frqs
> +    etaneg = etaneg_part / cpmg_frqs
> +
> +    # The mD value.
> +    mD = isqrt_pApBkex2 / (dpos * zpos) * (zpos + 
> 2.0*dw*sin(zpos*tcp)/sin((dpos + zpos)*tcp))
> +
> +    # The mZ value.
> +    mZ = -isqrt_pApBkex2 / (dneg * zneg) * (dneg - 
> 2.0*dw*sin(dneg*tcp)/sin((dneg + zneg)*tcp))
> +
> +    # The Q value.
> +    Q = 1 - mD**2 + mD*mZ - mZ**2 + 0.5*(mD + mZ)*sqrt_pBpA
> +    Q = Q.real
> +
> +    # The first eigenvalue.
> +    lambda1 = r20_kex - cpmg_frqs * arccosh(Dpos * cosh(etapos) - Dneg * 
> cos(etaneg))
> +
> +    # The full formula.
> +    R2eff = lambda1.real - inv_tcpmg * log(Q)
> +
> +    # 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)
> +
> +    # Parse back the value to update the back_calc class object.
>      for i in range(num_points):
> -        # Alias delta.
> -        delta = tcp[i]
> -
> -        # The full eta+/- values.
> -        etapos = etapos_part / cpmg_frqs[i]
> -        etaneg = etaneg_part / cpmg_frqs[i]
> -
> -        # The mD value.
> -        mD = isqrt_pApBkex2 / (dpos * zpos) * (zpos + 
> 2.0*dw*sin(zpos*delta)/sin((dpos + zpos)*delta))
> -
> -        # The mZ value.
> -        mZ = -isqrt_pApBkex2 / (dneg * zneg) * (dneg - 
> 2.0*dw*sin(dneg*delta)/sin((dneg + zneg)*delta))
> -
> -        # The Q value.
> -        Q = 1 - mD**2 + mD*mZ - mZ**2 + 0.5*(mD + mZ)*sqrt_pBpA
> -        Q = Q.real
> -
> -        # The first eigenvalue.
> -        lambda1 = r20_kex - cpmg_frqs[i] * arccosh(Dpos * cosh(etapos) - 
> Dneg * cos(etaneg))
> -
> -        # The full formula.
> -        back_calc[i] = lambda1.real - inv_tcpmg * log(Q)
> +        back_calc[i] = R2eff[i]
>
>
> _______________________________________________
> 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