Hi Troels, On top of the bug described at http://article.gmane.org/gmane.science.nmr.relax.devel/5938, there is a second bug here. I highly recommend reverting this commit.
The second bug is the deletion of the gamma < 0.0. This is absolutely essential, as later the square root of gamma is calculated. So if gamma is negative, which actually sometimes happens, then this code fails hard. That is the reason why I introduced this check. And because you have a square root of a negative number, which is not valid as these TAP03 equations are not designed for complex numbers, I decided that R2eff would be 1e100 to strongly penalise the optimisation when it reaches a region of the optimisation space where gamma < 0.0. Cheers, Edward On 26 May 2014 23:09, <[email protected]> wrote: > Author: tlinnet > Date: Mon May 26 23:09:52 2014 > New Revision: 23443 > > URL: http://svn.gna.org/viewcvs/relax?rev=23443&view=rev > Log: > Critical fix for the math domain catching of model TAP03. > > This was discovered with the added 8 unit tests demonstrating edge case 'no > Rex' failures. > > This follows from the ideas in the post > http://article.gmane.org/gmane.science.nmr.relax.devel/5858. > This is related to: task #7793: (https://gna.org/task/?7793) Speed-up of > dispersion models. > > This is to implement catching of math domain errors, before they occur. > > 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=23443&r1=23442&r2=23443&view=diff > ============================================================================== > --- branches/disp_speed/lib/dispersion/tap03.py (original) > +++ branches/disp_speed/lib/dispersion/tap03.py Mon May 26 23:09:52 2014 > @@ -60,7 +60,7 @@ > """ > > # Python module imports. > -from numpy import abs, arctan2, array, isfinite, min, sin, sqrt, sum > +from numpy import arctan2, array, 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, > num_points=None): > @@ -114,10 +114,6 @@ > > gamma = 1.0 + phi_ex*(sigma2 - kex2 + spin_lock_fields2) / (sigma2 + > kex2 + spin_lock_fields2)**2 > > - # Bad gamma. > - if min(gamma) < 0.0: > - return array([1e100]*num_points) > - > # 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. > @@ -135,17 +131,12 @@ > > # Catch zeros (to avoid pointless mathematical operations). > # This will result in no exchange, returning flat lines. > - if min(numer) == 0.0: > - return R1_cos_theta2 + R1rho_prime_sin_theta2 > + if numer == 0.0: > + return array([r1rho_prime]*num_points) > > # Denominator. > denom = waeff2*wbeff2/weff2 + kex2 - 2.0*hat_sin_theta2*phi_ex + (1.0 - > gamma)*spin_lock_fields2 > > - # Catch math domain error of dividing with 0. > - # This is when denom =0. > - if min(abs(denom)) == 0: > - return array([1e100]*num_points) > - > # R1rho calculation. > R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + hat_sin_theta2 * numer > / denom / gamma > > @@ -154,4 +145,4 @@ > if not isfinite(sum(R1rho)): > R1rho = array([1e100]*num_points) > > - return R1rho > + return R1rho > > > _______________________________________________ > 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
