Re: Error checking in the lib.dispersion modules (Re: r23270 - in /branches/disp_speed: lib/dispersion/tp02.py target_functions/relax_disp.py). -- May 21, 2014 - 08:49

Hi Troels,

Actually, if you add that check to the start of all lib.dispersion modules:

# Catch parameter values that will result in no exchange, returning
flat R2eff = R20 lines.
if dw == 0.0 or pA == 1.0 or kex == 0.0:
    return array([R20]*num_points)

it would then be worth checking if the other tests are necessary.
This one check might eliminate the need for most of the other checks
and the less checks there are, the faster the code will be.

Regards,

Edward



On 21 May 2014 08:37, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:
> Hi Troels,
>
> With some of these, you have to be very careful that the parameter
> combination actually does produce infinite R2eff values.  I think in
> this combination, it might not be the case.  What input parameter
> value causes weff2 to be zero?  To properly check, one needs to take
> the original equation and remove the parameter that is zero,
> collapsing the equations.  What looks like a divide by zero in the
> original equation might actually not be because of the equation
> simplifications.  For example in almost all cases of dw being zero,
> the R2eff equations should collapse to R20.  Or in this R1rho case, to
> R1.cos^2(theta) + R1rho'.sin^2(theta).  I.e. dw=0 means there is no
> exchange.  So you should never set R2eff to 1e100 in such a case, as
> you will no longer be able to produce the flat exchange-free R2eff
> lines.
>
> Note, you'll see that in some of my original checks in the
> lib.dispersion modules, when I am looking for a numerator that is zero
> before looking for a denominator that is zero.  You can see that in
> the original lib.dispersion.tp02 module.  But check has been removed
> from the disp_speed branch.  Without such checks, the dispersion
> models will incorrectly model no exchange as being almost infinite
> exchange.  It would be worth looking at all checks that I carefully
> added (they are still in the trunk), because every last one was very
> carefully thought out with the simplification it would make to the
> equations.  So every check should be in the disp_speed branch too.
> The checks I have added are not complete though.
>
> Maybe there is a far better way to perform some of these checks.  For
> example, catching parameter values that will result in zero exchange:
>
> # Catch parameter values that will result in no exchange, returning
> flat R2eff = R20 lines.
> if dw == 0.0 or pA == 1.0 or kex == 0.0:
>     return array([R20]*num_points)
>
> This could happen at the very start of all of the lib.dispersion
> modules.  Such a check would eliminate many of the numpy warnings
> where R2eff values end up being Inf, NaN or zero.  Actually, I feel
> rather stupid now for not having had thought of such a brilliant check
> before :S
>
> Regards,
>
> Edward
>
>
>
> On 20 May 2014 22:29,  <tlinnet@xxxxxxxxxxxxx> wrote:
> > Author: tlinnet
> > Date: Tue May 20 22:29:40 2014
> > New Revision: 23270
> >
> > URL: http://svn.gna.org/viewcvs/relax?rev=23270&view=rev
> > Log:
> > Math-domain catching for model TP02.
> >
> > task #7793: (https://gna.org/task/?7793) Speed-up of dispersion models.
> >
> > This is to implement catching of math domain errors, before they occur.
> > These can be found via the --numpy-raise function to the systemtests.
> >
> > To make the code look clean, the class object "back_calc" is no longer
> > being updated per time point, but is updated in the relax_disp target 
> > function in
> > one go.
> >
> > Modified:
> >     branches/disp_speed/lib/dispersion/tp02.py
> >     branches/disp_speed/target_functions/relax_disp.py
> >
> > Modified: branches/disp_speed/lib/dispersion/tp02.py
> > URL: 
> > http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/tp02.py?rev=23270&r1=23269&r2=23270&view=diff
> > ==============================================================================
> > --- branches/disp_speed/lib/dispersion/tp02.py  (original)
> > +++ branches/disp_speed/lib/dispersion/tp02.py  Tue May 20 22:29:40 2014
> > @@ -60,10 +60,10 @@
> >  """
> >
> >  # Python module imports.
> > -from numpy import arctan2, array, isfinite, sin, sum
> > +from numpy import abs, arctan2, array, isfinite, min, sin, sum
> >
> >
> > -def r1rho_TP02(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, back_calc=None, num_points=None):
> > +def r1rho_TP02(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):
> >      """Calculate the R1rho' values for the TP02 model.
> >
> >      See the module docstring for details.  This is the Trott and Palmer 
> > (2002) equation according to Korzhnev (J. Biomol. NMR (2003), 26, 39-48).
> > @@ -89,9 +89,7 @@
> >      @type spin_lock_fields:     numpy rank-1 float array
> >      @keyword spin_lock_fields2: The R1rho spin-lock field strengths 
> > squared (in rad^2.s^-2).  This is for speed.
> >      @type spin_lock_fields2:    numpy rank-1 float array
> > -    @keyword back_calc:         The array for holding the back calculated 
> > R1rho values.  Each element corresponds to one of the spin-lock fields.
> > -    @type back_calc:            numpy rank-1 float array
> > -    @keyword num_points:        The number of points on the dispersion 
> > curve, equal to the length of the spin_lock_fields and back_calc arguments.
> > +    @keyword num_points:        The number of points on the dispersion 
> > curve, equal to the length of the spin_lock_fields.
> >      @type num_points:           int
> >      """
> >
> > @@ -115,6 +113,13 @@
> >      wbeff2 = spin_lock_fields2 + db2       # Effective field at B.
> >      weff2 = spin_lock_fields2 + d2         # Effective field at 
> > pop-average.
> >
> > +    # Catch math domain error of dividing with 0.
> > +    # This is when weff2 = 0.
> > +    if min(abs(weff2)) == 0:
> > +        R2eff = array([1e100]*num_points)
> > +
> > +        return R2eff
> > +
> >      # The rotating frame flip angle.
> >      theta = arctan2(spin_lock_fields, d)
> >
> > @@ -127,6 +132,13 @@
> >      denom = waeff2 * wbeff2 / weff2 + kex2
> >      #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2
> >
> > +    # Catch math domain error of dividing with 0.
> > +    # This is when denom=0.
> > +    if min(abs(denom)) == 0:
> > +        R1rho = array([1e100]*num_points)
> > +
> > +        return R1rho
> > +
> >      # R1rho calculation.
> >      R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer / 
> > denom
> >
> > @@ -135,6 +147,4 @@
> >      if not isfinite(sum(R1rho)):
> >          R1rho = array([1e100]*num_points)
> >
> > -    # Parse back the value to update the back_calc class object.
> > -    for i in range(num_points):
> > -        back_calc[i] = R1rho[i]
> > +    return R1rho
> >
> > Modified: branches/disp_speed/target_functions/relax_disp.py
> > URL: 
> > http://svn.gna.org/viewcvs/relax/branches/disp_speed/target_functions/relax_disp.py?rev=23270&r1=23269&r2=23270&view=diff
> > ==============================================================================
> > --- branches/disp_speed/target_functions/relax_disp.py  (original)
> > +++ branches/disp_speed/target_functions/relax_disp.py  Tue May 20 
> > 22:29:40 2014
> > @@ -1875,7 +1875,7 @@
> >                  # Loop over the offsets.
> >                  for oi in range(self.num_offsets[0][si][mi]):
> >                      # Back calculate the R1rho values.
> > -                    r1rho_TP02(r1rho_prime=R20[r20_index], 
> > omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], 
> > pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi], 
> > spin_lock_fields=self.spin_lock_omega1[0][mi][oi], 
> > spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], 
> > back_calc=self.back_calc[0][si][mi][oi], 
> > num_points=self.num_disp_points[0][si][mi][oi])
> > +                    self.back_calc[0][si][mi][oi] = 
> > r1rho_TP02(r1rho_prime=R20[r20_index], 
> > omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], 
> > pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi], 
> > spin_lock_fields=self.spin_lock_omega1[0][mi][oi], 
> > spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], 
> > num_points=self.num_disp_points[0][si][mi][oi])
> >
> >                      # For all missing data points, set the 
> > back-calculated value to the measured values so that it has no effect on 
> > the chi-squared value.
> >                      for di in range(self.num_disp_points[0][si][mi][oi]):
> >
> >
> > _______________________________________________
> > 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