Hi Troels, Again for this dispersion function, you can use the zero 'denom' check to avoid all numpy warnings (http://thread.gmane.org/gmane.science.nmr.relax.scm/20966).
Regards, Edward On 19 May 2014 03:20, <[email protected]> wrote: > Author: tlinnet > Date: Mon May 19 03:20:43 2014 > New Revision: 23222 > > URL: http://svn.gna.org/viewcvs/relax?rev=23222&view=rev > Log: > Huge speed-up for model TAP03. > > task #7793: (https://gna.org/task/?7793) Speed-up of dispersion models. > > The change for running systemtest is: > test_tp02_data_to_tap03 > 13.869s -> 7.263s > > This is won by not checking single values in the R1rho array for math domain > errors, but calculating all steps, and in one single round check for finite > values. > If just one non-finite value is found, the whole array is returned with a > large > penalty of 1e100. > > This makes all calculations be the fastest numpy array way. > > Modified: > branches/disp_speed/lib/dispersion/tap03.py > > Modified: branches/disp_speed/lib/dispersion/tap03.py > URL: > http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/tap03.py?rev=23222&r1=23221&r2=23222&view=diff > ============================================================================== > --- branches/disp_speed/lib/dispersion/tap03.py (original) > +++ branches/disp_speed/lib/dispersion/tap03.py Mon May 19 03:20:43 2014 > @@ -60,7 +60,7 @@ > """ > > # Python module imports. > -from math import atan2, sin, sqrt > +from numpy import arctan2, isfinite, sin, sqrt, sum > > > def r1rho_TAP03(r1rho_prime=None, omega=None, offset=None, pA=None, pB=None, > dw=None, kex=None, R1=0.0, spin_lock_fields=None, spin_lock_fields2=None, > back_calc=None, num_points=None): > @@ -105,50 +105,42 @@ > phi_ex = pA * pB * dw**2 > numer = phi_ex * kex > > - # Loop over the dispersion points, back calculating the R1rho values. > + # The factors. > + da = Wa - offset # Offset of spin-lock from A. > + db = Wb - offset # Offset of spin-lock from B. > + d = W - offset # Offset of spin-lock from pop-average. > + > + # The gamma factor. > + sigma = pB*da + pA*db > + sigma2 = sigma**2 > + gamma = 1.0 + phi_ex*(sigma2 - kex2 + spin_lock_fields2) / (sigma2 + > kex2 + spin_lock_fields2)**2 > + > + # Special omega values. > + waeff2 = gamma*spin_lock_fields2 + da**2 # Effective field at A. > + wbeff2 = gamma*spin_lock_fields2 + db**2 # Effective field at B. > + weff2 = gamma*spin_lock_fields2 + d**2 # Effective field at > pop-average. > + > + # The rotating frame flip angle. > + theta = arctan2(spin_lock_fields, d) > + hat_theta = arctan2(sqrt(gamma)*spin_lock_fields, d) > + > + # Repetitive calculations (to speed up calculations). > + sin_theta2 = sin(theta)**2 > + hat_sin_theta2 = sin(hat_theta)**2 > + R1_cos_theta2 = R1 * (1.0 - sin_theta2) > + R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2 > + > + # Denominator. > + denom = waeff2*wbeff2/weff2 + kex2 - 2.0*hat_sin_theta2*phi_ex + (1.0 - > gamma)*spin_lock_fields2 > + > + # R1rho calculation. > + R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + hat_sin_theta2 * numer > / denom / gamma > + > + # 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(R1rho)): > + R1rho = array([1e100]*num_points) > + > + # Parse back the value to update the back_calc class object. > for i in range(num_points): > - # The factors. > - da = Wa - offset # Offset of spin-lock from A. > - db = Wb - offset # Offset of spin-lock from B. > - d = W - offset # Offset of spin-lock from pop-average. > - > - # The gamma factor. > - sigma = pB*da + pA*db > - sigma2 = sigma**2 > - gamma = 1.0 + phi_ex*(sigma2 - kex2 + spin_lock_fields2[i]) / > (sigma2 + kex2 + spin_lock_fields2[i])**2 > - > - # Bad gamma. > - if gamma < 0.0: > - back_calc[i] = 1e100 > - continue > - > - # Special omega values. > - waeff2 = gamma*spin_lock_fields2[i] + da**2 # Effective field at > A. > - wbeff2 = gamma*spin_lock_fields2[i] + db**2 # Effective field at > B. > - weff2 = gamma*spin_lock_fields2[i] + d**2 # Effective field at > pop-average. > - > - # The rotating frame flip angle. > - theta = atan2(spin_lock_fields[i], d) > - hat_theta = atan2(sqrt(gamma)*spin_lock_fields[i], d) > - > - # Repetitive calculations (to speed up calculations). > - sin_theta2 = sin(theta)**2 > - hat_sin_theta2 = sin(hat_theta)**2 > - R1_cos_theta2 = R1 * (1.0 - sin_theta2) > - R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2 > - > - # Catch zeros (to avoid pointless mathematical operations). > - if numer == 0.0: > - back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2 > - continue > - > - # Denominator. > - denom = waeff2*wbeff2/weff2 + kex2 - 2.0*hat_sin_theta2*phi_ex + > (1.0 - gamma)*spin_lock_fields2[i] > - > - # Avoid divide by zero. > - if denom == 0.0: > - back_calc[i] = 1e100 > - continue > - > - # R1rho calculation. > - back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2 + > hat_sin_theta2 * numer / denom / gamma > + back_calc[i] = R1rho[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
