Hi Troels, Have a look at this commit to the trunk:
http://www.mail-archive.com/[email protected]/msg20984.html Clearly kex = 1e5 is one of these edge cases where all R2eff must be set to R20 (this actually depends a little bit on dw too). But we probably should not catch this in the lib.dispersion modules. It is the model itself which should model this correctly - and I know that some cannot. If you use this test for the other models, you'll see for yourself. Regards, Edward On 22 May 2014 19:32, Edward d'Auvergne <[email protected]> wrote: > 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
