Re: r23451 - /branches/disp_speed/lib/dispersion/ns_cpmg_2site_3d.py -- May 27, 2014 - 18:18

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 <tlinnet@xxxxxxxxxxxxx> 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 <edward@xxxxxxxxxxxxx>:
> > 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