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 <[email protected]>:
> 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 <[email protected]> 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 <[email protected]>:
>>> 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, <[email protected]> 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
>>>> [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
