Author: tlinnet Date: Tue May 20 22:29:47 2014 New Revision: 23274 URL: http://svn.gna.org/viewcvs/relax?rev=23274&view=rev Log: Math-domain catching for model MP05.
task #7793: (https://gna.org/task/?7793) Speed-up of dispersion models. This is to implement catching of math domain errors, before they occur. These can be found via the --numpy-raise function to the systemtests. To make the code look clean, the class object "back_calc" is no longer being updated per time point, but is updated in the relax_disp target function in one go. Modified: branches/disp_speed/lib/dispersion/mp05.py branches/disp_speed/target_functions/relax_disp.py Modified: branches/disp_speed/lib/dispersion/mp05.py URL: http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/mp05.py?rev=23274&r1=23273&r2=23274&view=diff ============================================================================== --- branches/disp_speed/lib/dispersion/mp05.py (original) +++ branches/disp_speed/lib/dispersion/mp05.py Tue May 20 22:29:47 2014 @@ -60,10 +60,10 @@ """ # Python module imports. -from numpy import arctan2, array, isfinite, sin, sum +from numpy import abs, arctan2, array, isfinite, min, sin, sum -def r1rho_MP05(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): +def r1rho_MP05(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, num_points=None): """Calculate the R1rho' values for the TP02 model. See the module docstring for details. This is the Miloushev and Palmer (2005) equation according to Korzhnev (J. Biomol. NMR (2003), 26, 39-48). @@ -89,9 +89,7 @@ @type spin_lock_fields: numpy rank-1 float array @keyword spin_lock_fields2: The R1rho spin-lock field strengths squared (in rad^2.s^-2). This is for speed. @type spin_lock_fields2: numpy rank-1 float array - @keyword back_calc: The array for holding the back calculated R1rho values. Each element corresponds to one of the spin-lock fields. - @type back_calc: numpy rank-1 float array - @keyword num_points: The number of points on the dispersion curve, equal to the length of the spin_lock_fields and back_calc arguments. + @keyword num_points: The number of points on the dispersion curve, equal to the length of the spin_lock_fields. @type num_points: int """ @@ -113,6 +111,13 @@ wbeff2 = spin_lock_fields2 + db**2 # Effective field at B. weff2 = spin_lock_fields2 + d**2 # Effective field at pop-average. + # Catch math domain error of dividing with 0. + # This is when weff2 = 0. + if min(abs(weff2)) == 0: + R2eff = array([1e100]*num_points) + + return R2eff + # The rotating frame flip angle. theta = arctan2(spin_lock_fields, d) @@ -123,8 +128,24 @@ # Denominator. waeff2_wbeff2 = waeff2*wbeff2 - fact = 1.0 + 2.0*kex2*(pA*waeff2 + pB*wbeff2) / (waeff2_wbeff2 + weff2*kex2) + fact_denom = waeff2_wbeff2 + weff2*kex2 + + # Catch math domain error of dividing with 0. + # This is when fact_denom = 0. + if min(abs(fact_denom)) == 0: + R2eff = array([1e100]*num_points) + + return R2eff + + fact = 1.0 + 2.0*kex2*(pA*waeff2 + pB*wbeff2) / fact_denom denom = waeff2_wbeff2/weff2 + kex2 - sin_theta2*phi_ex*(fact) + + # Catch math domain error of dividing with 0. + # This is when denom=0. + if min(abs(denom)) == 0: + R1rho = array([1e100]*num_points) + + return R1rho # R1rho calculation. R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer / denom @@ -134,7 +155,5 @@ 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): - back_calc[i] = R1rho[i] + return R1rho Modified: branches/disp_speed/target_functions/relax_disp.py URL: http://svn.gna.org/viewcvs/relax/branches/disp_speed/target_functions/relax_disp.py?rev=23274&r1=23273&r2=23274&view=diff ============================================================================== --- branches/disp_speed/target_functions/relax_disp.py (original) +++ branches/disp_speed/target_functions/relax_disp.py Tue May 20 22:29:47 2014 @@ -1262,7 +1262,7 @@ # Loop over the offsets. for oi in range(self.num_offsets[0][si][mi]): # Back calculate the R1rho values. - r1rho_MP05(r1rho_prime=R20[r20_index], omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi], spin_lock_fields=self.spin_lock_omega1[0][mi][oi], spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], back_calc=self.back_calc[0][si][mi][oi], num_points=self.num_disp_points[0][si][mi][oi]) + self.back_calc[0][si][mi][oi] = r1rho_MP05(r1rho_prime=R20[r20_index], omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi], spin_lock_fields=self.spin_lock_omega1[0][mi][oi], spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], num_points=self.num_disp_points[0][si][mi][oi]) # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. for di in range(self.num_disp_points[0][si][mi][oi]): _______________________________________________ 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
