Btw.
Mx = fabs(Mint[1] / pA)
if Mx <= 0.0 or isNaN(Mx):
fabs, is float absolute value.
It can never be < 0.
Best
Troels
2014-05-27 18:18 GMT+02:00 Edward d'Auvergne <[email protected]>:
> Ah, the answer to that is simple - this model fails for high kex
> values. Therefore the test_ns_cpmg_2site_3D_no_rex8() unit test
> checking for high kex is meaningless. You are also pushing things far
> too far with a pA value of 0.99. This is far too extreme, so 0.9
> would be better or even 0.95. At 0.99 you are asking for computer
> truncation artifacts to appear and destroy the model - and differently
> on 32 vs. 64-bit systems. Note that there is a bug in the calculation
> of the dw_frq value in these unit tests:
>
> # Calculate spin Larmor frequencies in 2pi.
> frqs = sfrq * 2 * pi
>
> # Convert dw from ppm to rad/s.
> dw_frq = dw * frqs
>
> Here you have 'sfrq' in Hz units. But if you have a look at the
> 'frqs' argument to the relaxation dispersion target function class,
> you will see that it is in MHz, not Hz. This is to speed up the
> calculations by pre-removing 1e6, so that you don't have to divide by
> 1e6 in the target functions to convert out of ppm units. So you need
> to divide dw_frq by 1e6.
>
> To debug the unit tests, try printing out self.R2eff before the
> assertAlmostEqual() check, so that you can see a dispersion curve.
> Then change the kex value in test_ns_cpmg_2site_3D_no_rex8(). Try the
> following in order [1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10,
> ..., 1e40]. This can be done inside a loop over range(40) with
> printouts and turning the assertAlmostEqual() check off (but having
> something at the end to make the test fail so you can see a printout).
> You should then understand why this model does not reach R20 as kex
> increases. Strangely, you might see what I see - that the high kex
> case is ~2.1 and not 2.0. Hmmm, maybe there is a bug somewhere!
> You'll also see this model fail horribly due to truncation artifacts
> at around 1.0e+16.
>
> Regards,
>
> Edward
>
>
> On 27 May 2014 17:24, Troels Emtekær Linnet <[email protected]> wrote:
>> Hi Edward.
>>
>> This wont work.
>>
>> I need to replace all values in back_calc if just one test is violated.
>> I tried it with the unit tests, but it cannot get solved.
>>
>> Best
>> Troeles
>>
>> 2014-05-27 10:22 GMT+02:00 Edward d'Auvergne <[email protected]>:
>>> Hi Troels,
>>>
>>> I have a suggestion for this change too, specifically:
>>>
>>> @@ -123,6 +129,8 @@
>>> # The next lines calculate the R2eff using a two-point
>>> approximation, i.e. assuming that the decay is mono-exponential.
>>> Mx = fabs(Mint[1] / pA)
>>> if Mx <= 0.0 or isNaN(Mx):
>>> - back_calc[i] = 1e99
>>> + for i in range(num_points):
>>> + back_calc[i] = r20a
>>> + return
>>> else:
>>> back_calc[i]= -inv_tcpmg * log(Mx)
>>>
>>> Here you have introduced a second loop over the index i inside an
>>> index i loop. This is not the best idea, instead try:
>>>
>>> # The next lines calculate the R2eff using a two-point
>>> approximation, i.e. assuming that the decay is mono-exponential.
>>> Mx = fabs(Mint[1] / pA)
>>> if Mx <= 0.0 or isNaN(Mx):
>>> back_calc[i] = r20a
>>> else:
>>> back_calc[i]= -inv_tcpmg * log(Mx)
>>>
>>> There is no need to return at this point, just let the loop continue
>>> harmlessly. Also, one hint for the commit message would be to explain
>>> why you changed the 1e99 value to r20a. I.e. that when kex is very
>>> large, Mx will be zero, and although the log of zero is not defined,
>>> this is in the 'no Rex' region so it should return r20a rather than
>>> the -1/tcpmg*log(0.0) value of -1/tcpmg*-infinity (which is positive
>>> infinity). This makes me wonder if this numeric model is any good for
>>> large kex values which actually do cause exchange!
>>>
>>> Cheers,
>>>
>>> Edward
>>>
>>> On 27 May 2014 00:08, <[email protected]> wrote:
>>>> Author: tlinnet
>>>> Date: Tue May 27 00:08:08 2014
>>>> New Revision: 23451
>>>>
>>>> URL: http://svn.gna.org/viewcvs/relax?rev=23451&view=rev
>>>> Log:
>>>> Critical fix for the math domain catching of model NS CPMG 2site 3D.
>>>>
>>>> This was discovered with the added 8 unit tests demonstrating edge case
>>>> 'no Rex' failures.
>>>>
>>>> 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.
>>>>
>>>> Modified:
>>>> branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py
>>>>
>>>> Modified: branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py
>>>> URL:
>>>> http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py?rev=23451&r1=23450&r2=23451&view=diff
>>>> ==============================================================================
>>>> --- branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py (original)
>>>> +++ branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py Tue May 27
>>>> 00:08:08 2014
>>>> @@ -103,6 +103,12 @@
>>>> @type power: numpy int16, rank-1 array
>>>> """
>>>>
>>>> + # Catch parameter values that will result in no exchange, returning
>>>> flat R2eff = R20 lines (when kex = 0.0, k_AB = 0.0).
>>>> + if dw == 0.0 or pA == 1.0 or k_AB == 0.0:
>>>> + for i in range(num_points):
>>>> + back_calc[i] = r20a
>>>> + return
>>>> +
>>>> # The matrix R that contains all the contributions to the evolution,
>>>> i.e. relaxation, exchange and chemical shift evolution.
>>>> R = rcpmg_3d(R1A=r10a, R1B=r10b, R2A=r20a, R2B=r20b, pA=pA, pB=pB,
>>>> dw=dw, k_AB=k_AB, k_BA=k_BA)
>>>>
>>>> @@ -123,6 +129,8 @@
>>>> # The next lines calculate the R2eff using a two-point
>>>> approximation, i.e. assuming that the decay is mono-exponential.
>>>> Mx = fabs(Mint[1] / pA)
>>>> if Mx <= 0.0 or isNaN(Mx):
>>>> - back_calc[i] = 1e99
>>>> + for i in range(num_points):
>>>> + back_calc[i] = r20a
>>>> + return
>>>> else:
>>>> back_calc[i]= -inv_tcpmg * log(Mx)
>>>>
>>>>
>>>> _______________________________________________
>>>> 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
