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
