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, <tlinnet@xxxxxxxxxxxxx> 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
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
Re: r23451 - /branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py