Oh, that must be because the checks have been copied from one of the other numeric models. Deleting it should slightly increase the speed of the model. The fabs() part comes from the fitting_main_kex.py script from Mathilde Lescanne, Paul Schanda, and Dominique Marion (see http://thread.gmane.org/gmane.science.nmr.relax.devel/4138, https://gna.org/task/?7712#comment2 and https://gna.org/support/download.php?file_id=18262). So maybe they introduced it to avoid the negative values. I don't think the absolute value is from some original formula, so it is equivalent to my Mx <= 0.0 check, i.e. just added to avoid this problem (the result is quite different though).
Regards, Edward On 27 May 2014 19:21, Troels Emtekær Linnet <[email protected]> wrote: > Btw. > > Mx = fabs(Mint[1] / pA) > if Mx <= 0.0 or isNaN(Mx): > > fabs, is float absolute value. > It can never be < 0. > > Best > Troels > > > 2014-05-27 18:18 GMT+02:00 Edward d'Auvergne <[email protected]>: >> Ah, the answer to that is simple - this model fails for high kex >> values. Therefore the test_ns_cpmg_2site_3D_no_rex8() unit test >> checking for high kex is meaningless. You are also pushing things far >> too far with a pA value of 0.99. This is far too extreme, so 0.9 >> would be better or even 0.95. At 0.99 you are asking for computer >> truncation artifacts to appear and destroy the model - and differently >> on 32 vs. 64-bit systems. Note that there is a bug in the calculation >> of the dw_frq value in these unit tests: >> >> # Calculate spin Larmor frequencies in 2pi. >> frqs = sfrq * 2 * pi >> >> # Convert dw from ppm to rad/s. >> dw_frq = dw * frqs >> >> Here you have 'sfrq' in Hz units. But if you have a look at the >> 'frqs' argument to the relaxation dispersion target function class, >> you will see that it is in MHz, not Hz. This is to speed up the >> calculations by pre-removing 1e6, so that you don't have to divide by >> 1e6 in the target functions to convert out of ppm units. So you need >> to divide dw_frq by 1e6. >> >> To debug the unit tests, try printing out self.R2eff before the >> assertAlmostEqual() check, so that you can see a dispersion curve. >> Then change the kex value in test_ns_cpmg_2site_3D_no_rex8(). Try the >> following in order [1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, >> ..., 1e40]. This can be done inside a loop over range(40) with >> printouts and turning the assertAlmostEqual() check off (but having >> something at the end to make the test fail so you can see a printout). >> You should then understand why this model does not reach R20 as kex >> increases. Strangely, you might see what I see - that the high kex >> case is ~2.1 and not 2.0. Hmmm, maybe there is a bug somewhere! >> You'll also see this model fail horribly due to truncation artifacts >> at around 1.0e+16. >> >> Regards, >> >> Edward >> >> >> On 27 May 2014 17:24, Troels Emtekær Linnet <[email protected]> wrote: >>> Hi Edward. >>> >>> This wont work. >>> >>> I need to replace all values in back_calc if just one test is violated. >>> I tried it with the unit tests, but it cannot get solved. >>> >>> Best >>> Troeles >>> >>> 2014-05-27 10:22 GMT+02:00 Edward d'Auvergne <[email protected]>: >>>> 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, <[email protected]> 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 >>>>> [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 _______________________________________________ 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
