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 - 19:33

Hi Ed.

I think I have fixed the 1e100 to R20.
That was in commit 23297 - 23304

I still have to look at:
dpl94
lm63

And, I will look at the unit tests.

Best
Troels

2014-05-21 19:21 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
> Let me reinforce that again if you didn't quite get it yet ;)  In the
> 'disp_speed' branch, there are a lot of cases where an array of 1e100
> is returned when an array of R20 should be returned.  Unit tests would
> instantly show these.  If you put a catch at the start for all
> parameter values that result in no chemical exchange, then you will
> avoid most problems seen with the --numpy-raise flag, and the
> optimisation space you see with the dx.map user function might improve
> a lot.  And these single value checks will be faster than the numpy
> array checks further into the functions, most of which will become
> redundant.
>
> Regards,
>
> Edward
>
>
>
> On 21 May 2014 19:10, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:
> > Hi Troels,
> >
> > You should really consider adding this check to absolutely all
> > lib.dispersion modules.  And as close to the first line as possible
> > (for most it should be on line one, even if you have to replicate some
> > calculations just for the 'No Rex' case).  It should magically make
> > most problems go away!  Especially in the numeric dispersion models.
> > And then many checks will no longer be necessary
> > (https://gna.org/bugs/?22065 and most errors you see in the test suite
> > with --numpy-raise will disappear).
> >
> > Regards,
> >
> > Edward
> >
> >
> > On 21 May 2014 08:48, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:
> > > 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
>
> _______________________________________________
> 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