:) 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
