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, <[email protected]> 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
> [email protected]
>
> 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
[email protected]
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
