Re: r23220 - /branches/disp_speed/lib/dispersion/b14.py -- May 19, 2014 - 14:14

Hi Ed.

I can say, that I am in the process of re-inserting your code. :-)

Best
Troels

2014-05-19 14:09 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
> :)  Well, actually there is often only one failure point and I have
> identified those in most lib.dispersion.* modules for you already ;)
> I have spent a lot of time running these models, not only in the
> system tests but also in the test_suite/shared_data/dispersion/
> directories, finding exactly what these failure points are!  Therefore
> just look at the lib.dispersion.* modules in the relax trunk, and you
> will see the exact tests needed.  I spend a lot of time and effort
> identifying and creating these.  This hard work has been deleted in
> your disp_speed branch, but you can reintroduce it in a fast way.
> Most will be of the form "if x[i] == 0.0" which, for your branch, can
> be converted to "if numpy.min(numpy.abs(x)) == 0.0", replacing the "if
> not numpy.isfinite()" checks you have introduced.
>
> The isfinite() checks are too late in the calculation, hence you see a
> numpy warning.  The "== 0.0" check catches the problem before the
> numpy calculation which generates the values of infinity, so
> reintroducing these will solve the numpy problems that often confuse
> users.  I have already done the hard part for you :)
>
> Regards,
>
> Edward
>
>
>
> On 19 May 2014 13:15, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
> > Hi Ed.
> >
> > The problem is, that 2-4 lines of code in a lib file can trigger a
> > math domain error.
> >
> > To figure out what each function is "sad" about, one would need a
> > little inspection per lib file, and
> > the code would get full of checkings and handlings. That slows it down...
> >
> > It should be, that in 95 percent of the time, the minimise function is
> > just happy and calculating.
> > Only at boundaries values, the math domain errors can get triggered.
> >
> > Rather that 5% of the time, it calculates to much, that 95% is full of 
> > checking.
> >
> > Best
> > Troels
> >
> >
> >
> >
> > 2014-05-19 10:18 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
> > > Hi Troels,
> > >
> > > There is actually a better way that is just as fast.  That is to catch
> > > the value prior to the divide by zero, or what ever operation creates
> > > the Inf values.  Then you can fill back_calc with 1e100 and return.
> > > That way the test will pass when running with --numpy-raise, and the
> > > user will not see all the numpy warning messages.  Having the check
> > > earlier might even be a tiny bit faster.
> > >
> > > Regards,
> > >
> > > Edward
> > >
> > >
> > >
> > >
> > > On 19 May 2014 00:51,  <tlinnet@xxxxxxxxxxxxx> wrote:
> > > > Author: tlinnet
> > > > Date: Mon May 19 00:51:20 2014
> > > > New Revision: 23220
> > > >
> > > > URL: http://svn.gna.org/viewcvs/relax?rev=23220&view=rev
> > > > Log:
> > > > Huge speed-up of model B14.
> > > >
> > > > task #7793: (https://gna.org/task/?7793) Speed-up of dispersion models.
> > > >
> > > > Time test for systemtests:
> > > >
> > > > test_baldwin_synthetic
> > > > 2.626s -> 1.990s
> > > >
> > > > test_baldwin_synthetic_full
> > > > 18.326s -> 13.742s
> > > >
> > > > This is won by not checking single values in the R2eff array for math 
> > > > domain
> > > > errors, but calculating all steps, and in one single round check for 
> > > > finite values.
> > > > If just one non-finite value is found, the whole array is returned with 
> > > > a large
> > > > penalty of 1e100.
> > > >
> > > > This makes all calculations be the fastest numpy array way.
> > > >
> > > > Modified:
> > > >     branches/disp_speed/lib/dispersion/b14.py
> > > >
> > > > Modified: branches/disp_speed/lib/dispersion/b14.py
> > > > URL: 
> > > > http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/b14.py?rev=23220&r1=23219&r2=23220&view=diff
> > > > ==============================================================================
> > > > --- branches/disp_speed/lib/dispersion/b14.py   (original)
> > > > +++ branches/disp_speed/lib/dispersion/b14.py   Mon May 19 00:51:20 2014
> > > > @@ -110,8 +110,7 @@
> > > >  """
> > > >
> > > >  # Python module imports.
> > > > -import numpy
> > > > -from numpy import arccosh, array, cos, cosh, in1d, log, nonzero, sin, 
> > > > sinh, sqrt, power
> > > > +from numpy import arccosh, array, cos, cosh, isfinite, log, power, sin, 
> > > > sinh, sqrt, sum
> > > >
> > > >  # Repetitive calculations (to speed up calculations).
> > > >  g_fact = 1/sqrt(2)
> > > > @@ -216,70 +215,31 @@
> > > >      # Real. The v_1c in paper.
> > > >      v1c = F0 * cosh(E0) - F2 * cos(E2)
> > > >
> > > > -    # Catch devision with zero in y, when v3 = 0. v3 is 0, when v1c = 1.
> > > > -    # If no 1.0, perform normally.
> > > > -    if not in1d(1.0, v1c):
> > > > -        # Exact result for v2v3.
> > > > -        v3 = sqrt(v1c**2 - 1.)
> > > > -
> > > > -        y = power( (v1c - v3) / (v1c + v3), ncyc)
> > > > -
> > > > -        Tog = 0.5 * (1. + y) + (1. - y) * v5 / (2. * v3 * N )
> > > > -
> > > > -        # Find where Tog has negative values.
> > > > -        neg_index = nonzero(Tog.real < 0.0)[0]
> > > > -
> > > > -        # Do normal calculation
> > > > -        if len(neg_index) == 0:
> > > > -            ## -1/Trel * log(LpreDyn).
> > > > -            # Rpre = (r20a + r20b + kex) / 2.0
> > > > -
> > > > -            ## Carver and Richards (1972)
> > > > -            # R2eff_CR72 = Rpre - inv_tcpmg * ncyc *  arccosh(v1c.real)
> > > > -
> > > > -            ## Baldwin final.
> > > > -            # Estimate R2eff. relax_time = Trel = 1/inv_tcpmg.
> > > > -            # R2eff = R2eff_CR72 - inv_tcpmg * log(Tog.real)
> > > > -
> > > > -            # Fastest calculation.
> > > > -            R2eff = (r20a + r20b + kex) / 2.0  - inv_tcpmg * ( ncyc *  
> > > > arccosh(v1c.real) + log(Tog.real) )
> > > > -
> > > > -            # Loop over the time points, back calculating the R2eff 
> > > > values.
> > > > -            for i in range(num_points):
> > > > -
> > > > -                # Put values back.
> > > > -                back_calc[i] = R2eff[i]
> > > > -
> > > > -        else:
> > > > -            # Loop over each point.
> > > > -            for i in range(num_points):
> > > > -
> > > > -                # Return large value
> > > > -                if i in neg_index:
> > > > -                    back_calc[i] = 1e100
> > > > -
> > > > -                else:
> > > > -                    v3 = sqrt(v1c[i]**2 - 1.)
> > > > -                    y = power( (v1c[i] - v3) / (v1c[i] + v3), ncyc[i])
> > > > -                    Tog = 0.5 * (1. + y) + (1. - y) * v5[i] / (2. * v3 
> > > > * N )
> > > > -                    R2eff = (r20a + r20b + kex) / 2.0  - inv_tcpmg * ( 
> > > > ncyc[i] *  arccosh(v1c[i].real) + log(Tog.real) )
> > > > -                    back_calc[i] = R2eff
> > > > -
> > > > -    # This section is for catching math domain errors.
> > > > -    else:
> > > > -        # Find index where
> > > > -        one_indexes = nonzero(v1c == 1.0)[0]
> > > > -
> > > > -        # Loop over each point.
> > > > -        for i in range(num_points):
> > > > -
> > > > -            # Return large value
> > > > -            if i in one_indexes:
> > > > -                back_calc[i] = 1e100
> > > > -
> > > > -            else:
> > > > -                v3 = sqrt(v1c[i]**2 - 1.)
> > > > -                y = power( (v1c[i] - v3) / (v1c[i] + v3), ncyc[i])
> > > > -                Tog = 0.5 * (1. + y) + (1. - y) * v5[i] / (2. * v3 * N )
> > > > -                R2eff = (r20a + r20b + kex) / 2.0  - inv_tcpmg * ( 
> > > > ncyc[i] *  arccosh(v1c[i].real) + log(Tog.real) )
> > > > -                back_calc[i] = R2eff
> > > > +    # Exact result for v2v3.
> > > > +    v3 = sqrt(v1c**2 - 1.)
> > > > +
> > > > +    y = power( (v1c - v3) / (v1c + v3), ncyc)
> > > > +
> > > > +    Tog = 0.5 * (1. + y) + (1. - y) * v5 / (2. * v3 * N )
> > > > +
> > > > +    ## -1/Trel * log(LpreDyn).
> > > > +    # Rpre = (r20a + r20b + kex) / 2.0
> > > > +
> > > > +    ## Carver and Richards (1972)
> > > > +    # R2eff_CR72 = Rpre - inv_tcpmg * ncyc *  arccosh(v1c.real)
> > > > +
> > > > +    ## Baldwin final.
> > > > +    # Estimate R2eff. relax_time = Trel = 1/inv_tcpmg.
> > > > +    # R2eff = R2eff_CR72 - inv_tcpmg * log(Tog.real)
> > > > +
> > > > +    # Fastest calculation.
> > > > +    R2eff = (r20a + r20b + kex) / 2.0  - inv_tcpmg * ( ncyc *  
> > > > arccosh(v1c.real) + log(Tog.real) )
> > > > +
> > > > +    # Catch errors, taking a sum over array is the fastest way to check 
> > > > for
> > > > +    # +/- inf (infinity) and nan (not a number).
> > > > +    if not isfinite(sum(R2eff)):
> > > > +        R2eff = array([1e100]*num_points)
> > > > +
> > > > +    # Parse back the value to update the back_calc class object.
> > > > +    for i in range(num_points):
> > > > +        back_calc[i] = R2eff[i]
> > > >
> > > >
> > > > _______________________________________________
> > > > 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