Anyway, the test_ns_cpmg_2site_3D_no_rex8() test should be deleted as
the model does not function correctly for large kex values. You have
forced the correct behaviour when setting kex to 1e20, but it is not
correct. If I print out all Mx values for this kex value for this
test (deactivating the Mx <= 0.0 check), I see:
5.63779126881e+24
-2.13549768603e+23
1.86632678256e+27
485.247193505
485.247193505
485.247193505
-1.76141650845e+196
This is a total numerical collapse of the model. The fact that two of
these numbers being negative means that the check:
if Mx <= 0.0 or isNaN(Mx):
for i in range(num_points):
back_calc[i] = r20a
return
forces this failure point to be r20a. So the nested double loop "for
i in range(num_points)" is hiding a problem, not fixing a problem.
This is not good. It might be possible to come up with a kex value
close to 1e20 where no Mx value is negative. The model has failed on
an extreme edge case, i.e. extremely fast motions - one millionth of a
femtosecond, and there is nothing you can do about that. This is how
the magnetisation evolution matrix has been constructed, and it is a
consequence of the matrix_exponential() function. We have to just
accept this behaviour and let the model fail hard, as it should! The
nested looping can also be deleted, as it is only required to allow
the test_ns_cpmg_2site_3D_no_rex8() unit test to pass.
Regards,
Edward
On 27 May 2014 20:55, Edward d'Auvergne <[email protected]> wrote:
> Oh, that must be because the checks have been copied from one of the
> other numeric models. Deleting it should slightly increase the speed
> of the model. The fabs() part comes from the fitting_main_kex.py
> script from Mathilde Lescanne, Paul Schanda, and Dominique Marion (see
> http://thread.gmane.org/gmane.science.nmr.relax.devel/4138,
> https://gna.org/task/?7712#comment2 and
> https://gna.org/support/download.php?file_id=18262). So maybe they
> introduced it to avoid the negative values. I don't think the
> absolute value is from some original formula, so it is equivalent to
> my Mx <= 0.0 check, i.e. just added to avoid this problem (the result
> is quite different though).
>
> Regards,
>
> Edward
>
>
>
> On 27 May 2014 19:21, Troels Emtekær Linnet <[email protected]> wrote:
>> 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
