Hi Edward. Thank you for bringing in a very good demonstration data set. I do though see some problems in the setup. For spin 2, the error of the 800 MHz is way way to big. And also the graph "r1rho_vs_theta_inter_disp_2_N.png", also tells you, that you hit smack on chemical shift of residue 2, and the theta angle is not really covered. What generated the data: pA = 0.7654321 kex = 1e5 delta_omega = [7.0, 9.0] R2 (spin 1) = [10.0, 15.0] R2 (spin 2) = [12.0, 18.0] I tried a run, with the old linear constraints. ------------ With current setup Optimised parameter values: r2 (R1rho - 500.00000000 MHz) 9.999998314587026 r2 (R1rho - 800.00000000 MHz) 14.999995576824961 dw 7.000000871710442 pA 0.765432101642959 kex 99999.999999762905645 Optimised parameter values: r2 (R1rho - 500.00000000 MHz) 14.615528895956679 r2 (R1rho - 800.00000000 MHz) 24.381304162372505 dw 19.106783529541094 pA 0.966807181504147 kex 99999.999999999970896 chi2: 0.13578443834461 The spin cluster [':1@N']. # Data pipe Num_params_(k) Num_data_sets_(n) Chi2 Criterion No Rex - relax_disp 2 22 16229.31218 16233.31218 No Rex R1rho off res - relax_disp 2 22 221.50932 225.50932 TP02 - relax_disp 5 22 0.00000 10.00000 The model from the data pipe 'TP02 - relax_disp' has been selected. The spin cluster [':2@N']. # Data pipe Num_params_(k) Num_data_sets_(n) Chi2 Criterion No Rex - relax_disp 2 22 1.35267 5.35267 No Rex R1rho off res - relax_disp 2 22 1.45476 5.45476 TP02 - relax_disp 5 22 0.13578 10.13578 --------------- With old linear constraints What generated the data: pA = 0.7654321 kex = 1e5 delta_omega = [7.0, 9.0] R2 (spin 1) = [10.0, 15.0] R2 (spin 2) = [12.0, 18.0] Optimised parameter values: r2 (R1rho - 500.00000000 MHz) 9.999975946549126 r2 (R1rho - 800.00000000 MHz) 14.999938526871542 dw 7.000022273965653 pA 0.765432830477209 kex 100000.145859482014203 Optimised parameter values: r2 (R1rho - 500.00000000 MHz) 4.939869976036660 r2 (R1rho - 800.00000000 MHz) 0.000000000000000 dw 15.009280996951436 pA 0.861874627893427 kex 124724.354893419513246 chi2 : 0.000187103787660835 The spin cluster [':1@N']. # Data pipe Num_params_(k) Num_data_sets_(n) Chi2 Criterion No Rex - relax_disp 2 22 16229.31218 16233.31218 No Rex R1rho off res - relax_disp 2 22 221.50932 225.50932 TP02 - relax_disp 5 22 0.00000 10.00000 The model from the data pipe 'TP02 - relax_disp' has been selected. The spin cluster [':2@N']. # Data pipe Num_params_(k) Num_data_sets_(n) Chi2 Criterion No Rex - relax_disp 2 22 1.35267 5.35267 No Rex R1rho off res - relax_disp 2 22 1.45476 5.45476 TP02 - relax_disp 5 22 0.00019 10.00019 The model from the data pipe 'No Rex - relax_disp' has been selected. So both methods have problems. And any inspection of dataset for residue 2, would also rule it out. I am not sure, what you have tried to proof here. 2014-08-19 15:21 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi, 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