Hmmm, interesting. I might need some time thinking about the opposite type of unit test - the edge cases where R2eff should be Inf rather than R20. It should be easy enough to create a unit test with kex set to something big, say 1e6, and then see what happens with and without the check. But at some point you'll be fast enough that you won't see exchange any more, so really high kex values should also cause R2eff to be set to R20. Take an extreme example of tex = 1 ns.rad^-1 (i.e. kex = 1/tex = 1e8). This is the same time scale as the fast Lipari and Szabo model-free dynamics. So you would see something in a model-free analysis but absolutely nothing in a dispersion analysis. I wonder if our dispersion models capture this condition correctly?
Regards, Edward On 22 May 2014 19:19, Troels Emtekær Linnet <[email protected]> wrote: > I am more afraid for the cases when kex is 10000, or something similar. > > This was what happened in: > https://gna.org/bugs/index.php?22032 > bug #22032: Minimisation explodes when using Grid_INC=None > > Best > Troels > > > 2014-05-22 19:03 GMT+02:00 Edward d'Auvergne <[email protected]>: >> Hi Troels, >> >> Now that the 'no Rex' unit tests for this LM63 model have been >> implemented, I tried out the following in your branch: >> >> """ >> Index: lib/dispersion/lm63.py >> =================================================================== >> --- lib/dispersion/lm63.py (revision 23355) >> +++ lib/dispersion/lm63.py (working copy) >> @@ -99,11 +99,4 @@ >> kex_4 = 4.0 / kex >> >> # Calculate R2eff. >> - R2eff = r20 + rex * (1.0 - kex_4 * cpmg_frqs * tanh(kex / (4.0 * >> cpmg_frqs))) >> - >> - # Catch errors, taking a sum over array is the fastest way to check for >> - # +/- inf (infinity) and nan (not a number). >> - if not isfinite(sum(R2eff)): >> - R2eff = array([1e100]*num_points) >> - >> - return R2eff >> \ No newline at end of file >> + return r20 + rex * (1.0 - kex_4 * cpmg_frqs * tanh(kex / (4.0 * >> cpmg_frqs))) >> """ >> >> This change does not effect any system or unit test in the test suite, >> even with the --numpy-raise flag! The infinite condition appears to >> be gone now that there are 'no Rex' checks on the first lines of the >> R2eff calculating function. You can do your own testing and see. The >> catching of 'no Rex' at the start might allow us to remove many other >> checks, speeding up the code a little. >> >> Regards, >> >> Edward >> >> >> On 22 May 2014 18:27, <[email protected]> wrote: >>> Author: tlinnet >>> Date: Thu May 22 18:27:42 2014 >>> New Revision: 23349 >>> >>> URL: http://svn.gna.org/viewcvs/relax?rev=23349&view=rev >>> Log: >>> Added 7 unit tests demonstrating edge case 'no Rex' failures of the model >>> 'LM63'. >>> >>> 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. >>> >>> These tests cover all parameter value combinations which result in no >>> exchange: >>> >>> - dw = 0.0, >>> - pA = 1.0, >>> - kex = 0.0, >>> - dw = 0.0 and pA = 1.0, >>> - dw = 0.0 and kex = 0.0, >>> - pA = 1.0 and kex = 0.0, >>> - dw = 0.0, pA = 1.0, and kex = 0.0. >>> >>> Modified: >>> branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/__init__.py >>> branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/test_lm63.py >>> >>> Modified: >>> branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/__init__.py >>> URL: >>> http://svn.gna.org/viewcvs/relax/branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/__init__.py?rev=23349&r1=23348&r2=23349&view=diff >>> ============================================================================== >>> --- branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/__init__.py >>> (original) >>> +++ branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/__init__.py >>> Thu May 22 18:27:42 2014 >>> @@ -1,6 +1,7 @@ >>> >>> ############################################################################### >>> # >>> # >>> # Copyright (C) 2013-2014 Edward d'Auvergne >>> # >>> +# Copyright (C) 2013-2014 Troels E. Linnet >>> # >>> # >>> # >>> # This file is part of the program relax (http://www.nmr-relax.com). >>> # >>> # >>> # >>> @@ -23,5 +24,6 @@ >>> __all__ = [ >>> 'test___init__', >>> 'test_dpl94', >>> + 'test_lm63', >>> 'test_ns_cpmg_2site_expanded' >>> ] >>> >>> Modified: >>> branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/test_lm63.py >>> URL: >>> http://svn.gna.org/viewcvs/relax/branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/test_lm63.py?rev=23349&r1=23348&r2=23349&view=diff >>> ============================================================================== >>> --- branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/test_lm63.py >>> (original) >>> +++ branches/disp_speed/test_suite/unit_tests/_lib/_dispersion/test_lm63.py >>> Thu May 22 18:27:42 2014 >>> @@ -1,6 +1,7 @@ >>> >>> ############################################################################### >>> # >>> # >>> # Copyright (C) 2014 Edward d'Auvergne >>> # >>> +# Copyright (C) 2014 Troels E. Linnet >>> # >>> # >>> # >>> # This file is part of the program relax (http://www.nmr-relax.com). >>> # >>> # >>> # >>> @@ -20,15 +21,15 @@ >>> >>> ############################################################################### >>> >>> # Python module imports. >>> -from numpy import array, float64, int16, zeros >>> +from numpy import array, float64, int16, pi, zeros >>> from unittest import TestCase >>> >>> # relax module imports. >>> -from lib.dispersion.ns_cpmg_2site_expanded import >>> r2eff_ns_cpmg_2site_expanded >>> +from lib.dispersion.lm63 import r2eff_LM63 >>> >>> >>> -class Test_ns_cpmg_2site_expanded(TestCase): >>> - """Unit tests for the lib.dispersion.ns_cpmg_2site_expanded relax >>> module.""" >>> +class Test_lm63(TestCase): >>> + """Unit tests for the lib.dispersion.lm63 relax module.""" >>> >>> def setUp(self): >>> """Set up for all unit tests.""" >>> @@ -39,51 +40,59 @@ >>> self.dw = 0.5 >>> self.kex = 100.0 >>> >>> + # The spin Larmor frequencies. >>> + self.sfrq = 599.8908617*1E6 >>> + >>> # Required data structures. >>> self.num_points = 3 >>> - self.tcp = array([0.1, 0.2, 0.3], float64) >>> + self.cpmg_frqs = array([[2.5, 1.25, 0.83]], float64) >>> self.R2eff = zeros(3, float64) >>> - self.num_cpmg = array([1, 2, 3], int16) >>> >>> >>> def calc_r2eff(self): >>> """Calculate and check the R2eff values.""" >>> >>> # Parameter conversions. >>> - k_AB, k_BA = self.param_conversion(pA=self.pA, kex=self.kex) >>> + phi_ex_scaled = self.param_conversion(pA=self.pA, dw=self.dw, >>> sfrq=self.sfrq) >>> >>> # Calculate the R2eff values. >>> - R2eff = r2eff_ns_cpmg_2site_expanded(r20=self.r20, pA=self.pA, >>> dw=self.dw, k_AB=k_AB, k_BA=k_BA, relax_time=0.3, inv_relax_time=1/0.3, >>> tcp=self.tcp, num_points=self.num_points, num_cpmg=self.num_cpmg) >>> - >>> + R2eff = r2eff_LM63(r20=self.r20, phi_ex=phi_ex_scaled, >>> kex=self.kex, cpmg_frqs=self.cpmg_frqs, num_points=self.num_points) >>> # Check all R2eff values. >>> for i in range(self.num_points): >>> self.assertAlmostEqual(R2eff[i], self.r20) >>> >>> >>> - def param_conversion(self, pA=None, kex=None): >>> + def param_conversion(self, pA=None, dw=None, sfrq=None): >>> """Convert the parameters. >>> >>> @keyword pA: The population of state A. >>> @type pA: float >>> - @keyword kex: The rate of exchange. >>> - @type kex: float >>> - @return: The parameters {k_AB, k_BA}. >>> - @rtype: tuple of float >>> + @keyword dw: The chemical exchange difference between states A >>> and B in ppm. >>> + @type dw: float >>> + @keyword sfrq: The spin Larmor frequencies in Hz. >>> + @type sfrq: float >>> + @return: The parameters phi_ex_scaled >>> + @rtype: float >>> """ >>> >>> # Calculate pB. >>> pB = 1.0 - pA >>> >>> - # Exchange rates. >>> - k_BA = pA * kex >>> - k_AB = pB * kex >>> + # Calculate spin Larmor frequencies in 2pi. >>> + frqs = sfrq * 2 * pi >>> + >>> + # The phi_ex parameter value (pA * pB * delta_omega^2). >>> + phi_ex = pA * pB * dw**2 >>> + >>> + # Convert phi_ex from ppm^2 to (rad/s)^2. >>> + phi_ex_scaled = phi_ex * frqs**2 >>> >>> # Return all values. >>> - return k_AB, k_BA >>> + return phi_ex_scaled >>> >>> >>> - def test_ns_cpmg_2site_expanded_no_rex1(self): >>> - """Test the r2eff_ns_cpmg_2site_expanded() function for no >>> exchange when dw = 0.0.""" >>> + def test_lm63_no_rex1(self): >>> + """Test the r2eff_lm63() function for no exchange when dw = 0.0.""" >>> >>> # Parameter reset. >>> self.dw = 0.0 >>> @@ -92,8 +101,8 @@ >>> self.calc_r2eff() >>> >>> >>> - def test_ns_cpmg_2site_expanded_no_rex2(self): >>> - """Test the r2eff_ns_cpmg_2site_expanded() function for no >>> exchange when pA = 1.0.""" >>> + def test_lm63_no_rex2(self): >>> + """Test the r2eff_lm63() function for no exchange when pA = 1.0.""" >>> >>> # Parameter reset. >>> self.pA = 1.0 >>> @@ -102,8 +111,8 @@ >>> self.calc_r2eff() >>> >>> >>> - def test_ns_cpmg_2site_expanded_no_rex3(self): >>> - """Test the r2eff_ns_cpmg_2site_expanded() function for no >>> exchange when kex = 0.0.""" >>> + def test_lm63_no_rex3(self): >>> + """Test the r2eff_lm63() function for no exchange when kex = >>> 0.0.""" >>> >>> # Parameter reset. >>> self.kex = 0.0 >>> @@ -112,8 +121,8 @@ >>> self.calc_r2eff() >>> >>> >>> - def test_ns_cpmg_2site_expanded_no_rex4(self): >>> - """Test the r2eff_ns_cpmg_2site_expanded() function for no >>> exchange when dw = 0.0 and pA = 1.0.""" >>> + def test_lm63_no_rex4(self): >>> + """Test the r2eff_lm63() function for no exchange when dw = 0.0 >>> and pA = 1.0.""" >>> >>> # Parameter reset. >>> self.pA = 1.0 >>> @@ -123,8 +132,8 @@ >>> self.calc_r2eff() >>> >>> >>> - def test_ns_cpmg_2site_expanded_no_rex5(self): >>> - """Test the r2eff_ns_cpmg_2site_expanded() function for no >>> exchange when dw = 0.0 and kex = 0.0.""" >>> + def test_lm63_no_rex5(self): >>> + """Test the r2eff_lm63() function for no exchange when dw = 0.0 >>> and kex = 0.0.""" >>> >>> # Parameter reset. >>> self.dw = 0.0 >>> @@ -134,8 +143,8 @@ >>> self.calc_r2eff() >>> >>> >>> - def test_ns_cpmg_2site_expanded_no_rex6(self): >>> - """Test the r2eff_ns_cpmg_2site_expanded() function for no >>> exchange when pA = 1.0 and kex = 0.0.""" >>> + def test_lm63_no_rex6(self): >>> + """Test the r2eff_lm63() function for no exchange when pA = 1.0 >>> and kex = 0.0.""" >>> >>> # Parameter reset. >>> self.pA = 1.0 >>> @@ -145,8 +154,8 @@ >>> self.calc_r2eff() >>> >>> >>> - def test_ns_cpmg_2site_expanded_no_rex7(self): >>> - """Test the r2eff_ns_cpmg_2site_expanded() function for no >>> exchange when dw = 0.0, pA = 1.0, and kex = 0.0.""" >>> + def test_lm63_no_rex7(self): >>> + """Test the r2eff_lm63() function for no exchange when dw = 0.0, >>> pA = 1.0, and kex = 0.0.""" >>> >>> # Parameter reset. >>> self.dw = 0.0 >>> >>> >>> _______________________________________________ >>> 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
