Hi Troels, I have a suggestion for this change too, specifically: @@ -123,6 +129,8 @@ # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. Mx = fabs(Mint[1] / pA) if Mx <= 0.0 or isNaN(Mx): - back_calc[i] = 1e99 + for i in range(num_points): + back_calc[i] = r20a + return else: back_calc[i]= -inv_tcpmg * log(Mx) Here you have introduced a second loop over the index i inside an index i loop. This is not the best idea, instead try: # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. Mx = fabs(Mint[1] / pA) if Mx <= 0.0 or isNaN(Mx): back_calc[i] = r20a else: back_calc[i]= -inv_tcpmg * log(Mx) There is no need to return at this point, just let the loop continue harmlessly. Also, one hint for the commit message would be to explain why you changed the 1e99 value to r20a. I.e. that when kex is very large, Mx will be zero, and although the log of zero is not defined, this is in the 'no Rex' region so it should return r20a rather than the -1/tcpmg*log(0.0) value of -1/tcpmg*-infinity (which is positive infinity). This makes me wonder if this numeric model is any good for large kex values which actually do cause exchange! Cheers, Edward On 27 May 2014 00:08, <tlinnet@xxxxxxxxxxxxx> wrote:
Author: tlinnet Date: Tue May 27 00:08:08 2014 New Revision: 23451 URL: http://svn.gna.org/viewcvs/relax?rev=23451&view=rev Log: Critical fix for the math domain catching of model NS CPMG 2site 3D. 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/ns_cpmg_2site_3d.py Modified: branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py URL: http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py?rev=23451&r1=23450&r2=23451&view=diff ============================================================================== --- branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py (original) +++ branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py Tue May 27 00:08:08 2014 @@ -103,6 +103,12 @@ @type power: numpy int16, rank-1 array """ + # Catch parameter values that will result in no exchange, returning flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0). + if dw == 0.0 or pA == 1.0 or k_AB == 0.0: + for i in range(num_points): + back_calc[i] = r20a + return + # The matrix R that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution. R = rcpmg_3d(R1A=r10a, R1B=r10b, R2A=r20a, R2B=r20b, pA=pA, pB=pB, dw=dw, k_AB=k_AB, k_BA=k_BA) @@ -123,6 +129,8 @@ # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. Mx = fabs(Mint[1] / pA) if Mx <= 0.0 or isNaN(Mx): - back_calc[i] = 1e99 + for i in range(num_points): + back_calc[i] = r20a + return else: back_calc[i]= -inv_tcpmg * log(Mx) _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-commits mailing list relax-commits@xxxxxxx 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