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@gna.org 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