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 <[email protected]>: > :) 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 <[email protected]> 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 <[email protected]>: >>> 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, <[email protected]> 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 >>>> [email protected] >>>> >>>> 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 >>> [email protected] >>> >>> 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 _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list [email protected] 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
