Re: r23443 - /branches/disp_speed/lib/dispersion/tap03.py -- May 27, 2014 - 09:53

Hi Troels,

On top of the bug described at
http://article.gmane.org/gmane.science.nmr.relax.devel/5938, there is
a second bug here.  I highly recommend reverting this commit.

The second bug is the deletion of the gamma < 0.0.  This is absolutely
essential, as later the square root of gamma is calculated.  So if
gamma is negative, which actually sometimes happens, then this code
fails hard.  That is the reason why I introduced this check.  And
because you have a square root of a negative number, which is not
valid as these TAP03 equations are not designed for complex numbers, I
decided that R2eff would be 1e100 to strongly penalise the
optimisation when it reaches a region of the optimisation space where
gamma < 0.0.

Cheers,

Edward



On 26 May 2014 23:09,  <tlinnet@xxxxxxxxxxxxx> wrote:
> Author: tlinnet
> Date: Mon May 26 23:09:52 2014
> New Revision: 23443
>
> URL: http://svn.gna.org/viewcvs/relax?rev=23443&view=rev
> Log:
> Critical fix for the math domain catching of model TAP03.
>
> 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/tap03.py
>
> Modified: branches/disp_speed/lib/dispersion/tap03.py
> URL: 
> http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/tap03.py?rev=23443&r1=23442&r2=23443&view=diff
> ==============================================================================
> --- branches/disp_speed/lib/dispersion/tap03.py (original)
> +++ branches/disp_speed/lib/dispersion/tap03.py Mon May 26 23:09:52 2014
> @@ -60,7 +60,7 @@
>  """
>
>  # Python module imports.
> -from numpy import abs, arctan2, array, isfinite, min, sin, sqrt, sum
> +from numpy import arctan2, array, isfinite, sin, sqrt, sum
>
>
>  def r1rho_TAP03(r1rho_prime=None, omega=None, offset=None, pA=None, 
> pB=None, dw=None, kex=None, R1=0.0, spin_lock_fields=None, 
> spin_lock_fields2=None, num_points=None):
> @@ -114,10 +114,6 @@
>
>      gamma = 1.0 + phi_ex*(sigma2 - kex2 + spin_lock_fields2) / (sigma2 + 
> kex2 + spin_lock_fields2)**2
>
> -    # Bad gamma.
> -    if min(gamma) < 0.0:
> -        return array([1e100]*num_points)
> -
>      # Special omega values.
>      waeff2 = gamma*spin_lock_fields2 + da**2     # Effective field at A.
>      wbeff2 = gamma*spin_lock_fields2 + db**2     # Effective field at B.
> @@ -135,17 +131,12 @@
>
>      # Catch zeros (to avoid pointless mathematical operations).
>      # This will result in no exchange, returning flat lines.
> -    if min(numer) == 0.0:
> -        return R1_cos_theta2 + R1rho_prime_sin_theta2
> +    if numer == 0.0:
> +        return array([r1rho_prime]*num_points)
>
>      # Denominator.
>      denom = waeff2*wbeff2/weff2 + kex2 - 2.0*hat_sin_theta2*phi_ex + (1.0 
> - gamma)*spin_lock_fields2
>
> -    # Catch math domain error of dividing with 0.
> -    # This is when denom =0.
> -    if min(abs(denom)) == 0:
> -        return array([1e100]*num_points)
> -
>      # R1rho calculation.
>      R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + hat_sin_theta2 * 
> numer / denom / gamma
>
> @@ -154,4 +145,4 @@
>      if not isfinite(sum(R1rho)):
>          R1rho = array([1e100]*num_points)
>
> -    return R1rho
> +    return R1rho
>
>
> _______________________________________________
> 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