Anyway, the test_ns_cpmg_2site_3D_no_rex8() test should be deleted as the model does not function correctly for large kex values. You have forced the correct behaviour when setting kex to 1e20, but it is not correct. If I print out all Mx values for this kex value for this test (deactivating the Mx <= 0.0 check), I see: 5.63779126881e+24 -2.13549768603e+23 1.86632678256e+27 485.247193505 485.247193505 485.247193505 -1.76141650845e+196 This is a total numerical collapse of the model. The fact that two of these numbers being negative means that the check: if Mx <= 0.0 or isNaN(Mx): for i in range(num_points): back_calc[i] = r20a return forces this failure point to be r20a. So the nested double loop "for i in range(num_points)" is hiding a problem, not fixing a problem. This is not good. It might be possible to come up with a kex value close to 1e20 where no Mx value is negative. The model has failed on an extreme edge case, i.e. extremely fast motions - one millionth of a femtosecond, and there is nothing you can do about that. This is how the magnetisation evolution matrix has been constructed, and it is a consequence of the matrix_exponential() function. We have to just accept this behaviour and let the model fail hard, as it should! The nested looping can also be deleted, as it is only required to allow the test_ns_cpmg_2site_3D_no_rex8() unit test to pass. Regards, Edward On 27 May 2014 20:55, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:
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 <tlinnet@xxxxxxxxxxxxx> 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 <edward@xxxxxxxxxxxxx>: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 <tlinnet@xxxxxxxxxxxxx> 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 <edward@xxxxxxxxxxxxx>: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_______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list relax-devel@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-devel