Re: r25024 - /trunk/specific_analyses/relax_disp/parameters.pyRe: r25024 - /trunk/specific_analyses/relax_disp/parameters.pyHi, Please see r25074 (http://article.gmane.org/gmane.science.nmr.relax.scm/22824) for a demonstration of how the statement from this article is rubbish. Care must be taken when using statements in the literature, as there is a lot of conjecture published in papers which do not match reality. Many such statements also suffer as NMR advances and we can collect higher quality and greater quantities of data. The limits of the model are strongly influenced by the errors in the data, hence decreasing these errors extends the limits. If optimisation reaches such high kex values, a quick look at the dispersion curve will show that there is no dispersion. When kex is very, very high, this is identical to kex == 0.0. You could reset kex to 0.0 after the optimisation and re-perform the optimisation to see that it is almost equivalent. Therefore a value of say 20,000 tells you that there is no, or very little exchange (or that optimisation failed to find the minimum, though for the simplicity of the dispersion model space this is not really an issue). The Monte Carlo simulations and resultant kex errors will also tell you this. And finally model selection should choose the 'No Rex' model - again demonstrating that the high value is insignificant. If you put a constraint on kex to what you think a reasonable value is, then you are introducing bias into the optimisation. This is very bad - bias should be avoided at all cost. Here you create new fake new minimum by cutting the real minimum out of the space. But this new minimum has no relationship to physical reality. You must let the data tell you what the result is, and not the other way around. The upper constraint on kex is not to allow reasonable kex values to be found. It is to prevent an optimisation that would take forever to complete. When a parameter of the model to be optimised decides to shoot out to infinity - well the time for optimisation also decides to shoot out to infinity (though in most cases the optimisation algorithm dies before infinite time is reached). This is the one and only purpose for this constraint. The other purpose for constraints is to exclude regions of the optimisation space which are physically impossible, or are repetitive due to symmetries in the mathematical problem. Therefore I believe that this commit is not good practice. Regards, Edward On 19 August 2014 11:17, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
Hi Edward. Can you point to any reference, where kex is measured at higher values? I see no logic in having linear constraints for a search space higher than physically possible? The edge cases you mention (https://gna.org/bugs/?22024) is a very good example. I later learned for this case, that data was acquired at 3 different temperatures, to allow for robust interpretation of the data. Best Troels 2014-08-19 10:58 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:Hi Troels, I am unable to download that reference and the book with the reference in Goettingen is 500 km away from where I am now. Would you be able to send me a copy in a private mail? Cheers. I'm not convinced that this change should be made. The constraint of 1e6 was simply a large upper bound whereby no exchange should be visible. The constraint is used to prevent kex from heading out to infinity. So changing it to 1e4 or 1e5 will have no effect on this purpose of the constraint. What worries me here is that although exchange is very close to zero at these new limits, they are not exactly zero. This is especially the case for data with very, very small errors (test data with zero errors will show this). So the edge cases you have mapped out with OpenDX (https://gna.org/bugs/?22024 for example) will be highly effected by such a change. But this would be an edge case on the other side of the space. The new limits could chop a very weak minimum out of the optimisation space - i.e. you've put up a wall in front of these weak minima. This is a very weak edge case, but it nevertheless exists. As the purpose of the constraint is simply to stop kex from going to infinity, loosing these special edge case minima is not worth it. I would suggest reverting this change. Cheers, Edward On 15 August 2014 15:11, <tlinnet@xxxxxxxxxxxxx> wrote:Author: tlinnet Date: Fri Aug 15 15:11:00 2014 New Revision: 25024 URL: http://svn.gna.org/viewcvs/relax?rev=25024&view=rev Log: Modified the Linear Constraints for the exchange rates. For CPMG, the maximum kex should be 10^4, and for R1rho it should be 10^5. This is altered from the value of 10^6. The suggested restraints for 'kex' follows from article, on page 224: Nuclear Magnetic Resonance Methods for Quantifying Microsecond-to-Millisecond Motions in Biological Macromolecules. Palmer-III, Arthur G., Kroenke, Christopher D., Loria, J. Patrick Nucl. Magn. Reson. Biol. Macromol. B, 2001, Vol: 339, pages 204-238. U{DOI: 10.1016/S0076-6879(01)39315-1<http://dx.doi.org/10.1016/S0076-6879%2801%2939315-1>} Modified: trunk/specific_analyses/relax_disp/parameters.py Modified: trunk/specific_analyses/relax_disp/parameters.py URL: http://svn.gna.org/viewcvs/relax/trunk/specific_analyses/relax_disp/parameters.py?rev=25024&r1=25023&r2=25024&view=diff ============================================================================== --- trunk/specific_analyses/relax_disp/parameters.py (original) +++ trunk/specific_analyses/relax_disp/parameters.py Fri Aug 15 15:11:00 2014 @@ -34,7 +34,7 @@ from pipe_control import pipes from pipe_control.mol_res_spin import exists_mol_res_spin_data, return_spin from specific_analyses.relax_disp.data import count_spins, generate_r20_key, has_exponential_exp_type, loop_cluster, loop_exp_frq -from specific_analyses.relax_disp.variables import MODEL_LIST_MMQ, MODEL_M61B, MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR, MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR, PARAMS_R20 +from specific_analyses.relax_disp.variables import MODEL_LIST_ANALYTIC_R1RHO, MODEL_LIST_CPMG_ONLY, MODEL_LIST_MMQ, MODEL_LIST_NUMERIC_R1RHO, MODEL_M61B, MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR, MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR, PARAMS_R20 def assemble_param_vector(spins=None, key=None, sim_index=None): @@ -435,6 +435,12 @@ def linear_constraints(spins=None, scaling_matrix=None): """Set up the relaxation dispersion curve fitting linear constraint matrices A and b. + The suggested restraints for 'kex' follows from article, on page 224: + Nuclear Magnetic Resonance Methods for Quantifying Microsecond-to-Millisecond Motions in Biological Macromolecules. + Palmer-III, Arthur G., Kroenke, Christopher D., Loria, J. Patrick + Nucl. Magn. Reson. Biol. Macromol. B, 2001, Vol: 339, pages 204-238. + U{DOI: 10.1016/S0076-6879(01)39315-1<http://dx.doi.org/10.1016/S0076-6879%2801%2939315-1>}. + Standard notation ================= @@ -453,12 +459,12 @@ phi_ex_C >= 0 padw2 >= 0 dw >= 0 - 0 <= kex <= 2e6 - 0 <= k_AB <= 2e6 - 0 <= kB <= 2e6 - 0 <= kC <= 2e6 + 0 <= kex <= 1e4, for CPMG + 0 <= kex <= 1e5, for R1rho + 0 <= k_AB <= 1e4 + 0 <= kB <= 1e4 + 0 <= kC <= 1e4 tex >= 0 - k_AB >= 0 Matrix notation @@ -502,19 +508,22 @@ | | | | | | | 1 0 0 | | kex | | 0 | | | | | | | - |-1 0 0 | | kex | | -2e6 | + |-1 0 0 | | kex | |-1e4/-1e5| + | | | | | | + | 1 0 0 | | k_AB | | 0 | + | | | | | | + |-1 0 0 | | k_AB | | -1e4 | | | | | | | | 1 0 0 | | kB | | 0 | | | | | | | - |-1 0 0 | | kB | | -2e6 | + |-1 0 0 | | kB | | -1e4 | | | | | | | | 1 0 0 | | kC | | 0 | | | | | | | - |-1 0 0 | | kC | | -2e6 | + |-1 0 0 | | kC | | -1e4 | | | | | | | | 1 0 0 | | tex | | 0 | | | | | | | - | 1 0 0 | | k_AB | | 0 | @keyword spins: The list of spin data containers for the block. @@ -628,14 +637,21 @@ j += 1 break - # Exchange rates and times (0 <= k <= 2e6). + # Exchange rates and times (0 <= k <= 1e4) for CPMG and (0 <= k <= 1e5) for R1rho. elif param_name in ['kex', 'kex_AB', 'kex_AC', 'kex_BC', 'k_AB', 'kB', 'kC']: A.append(zero_array * 0.0) A.append(zero_array * 0.0) A[j][param_index] = 1.0 A[j+1][param_index] = -1.0 b.append(0.0) - b.append(-2e6 / scaling_matrix[param_index, param_index]) + # For CPMG experiments, (0 <= k <= 1e4). + if spins[0].model in MODEL_LIST_CPMG_ONLY + MODEL_LIST_MMQ: + b.append(-1e4 / scaling_matrix[param_index, param_index]) + # For R1rho experiments, (0 <= k <= 1e5). + elif spins[0].model in MODEL_LIST_ANALYTIC_R1RHO + MODEL_LIST_NUMERIC_R1RHO: + b.append(-1e5 / scaling_matrix[param_index, param_index]) + else: + b.append(-2e6 / scaling_matrix[param_index, param_index]) j += 2 # Exchange times (tex >= 0). _______________________________________________ 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