Very much so :) On a related note, you would get insane speed ups of the dispersion model if one or more of the loops in the target functions could be brought into these lib.dispersion modules and then this exact logic used. However I have never managed to make this work - though I haven't tried too hard. It would require quite large numpy data structures with multiply replicated data. Such a revolutionary idea would really need its own real subversion branch, and a lot of playing around testing ideas.
Regards, Edward On 8 May 2014 19:03, Troels Emtekær Linnet <tlin...@nmr-relax.com> wrote: > Yeah. > > I tried really hard to stay in numpy land as long as possible. > > Are you not here introducing the rather stupid mistake of always > calculating in a loop? > > I have seen the same for CR72. > > Best > Troels > > 2014-05-08 18:50 GMT+02:00 Edward d'Auvergne <edw...@nmr-relax.com>: >> Hi, >> >> For a test, and as a reference, I replaced all lines after the v5 = line >> with: >> >> """ >> # Real. The v_1c in paper. >> v1c = F0 * cosh(E0) - F2 * cos(E2) >> >> # Loop over each point. >> for i in range(num_points): >> # Catch devision with zero in y, when v3 = 0. v3 is 0, when v1c = 1. >> # If no 1.0, perform normally. >> if v1c[i] == 1: >> back_calc[i] = 1e100 >> continue >> >> # Exact result for v2v3. >> 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 ) >> >> # Find where Tog has negative values. >> if Tog.real < 0.0: >> back_calc[i] = 1e100 >> continue >> >> # Fastest calculation. >> back_calc[i] = (r20a + r20b + kex) / 2.0 - inv_tcpmg * ( >> ncyc[i] * arccosh(v1c[i].real) + log(Tog.real) ) >> """ >> >> This is far more compact and easier to read. However the >> Relax_disp.test_baldwin_synthetic_full system test time on my system >> went from 7.648s to 11.155s! This difference is because the old code >> uses the numpy speed advantage of operating on numpy arrays, and the >> above code does not. >> >> Regards, >> >> Edward >> >> On 8 May 2014 18:35, <tlin...@nmr-relax.com> wrote: >>> Author: tlinnet >>> Date: Thu May 8 18:35:51 2014 >>> New Revision: 23091 >>> >>> URL: http://svn.gna.org/viewcvs/relax?rev=23091&view=rev >>> Log: >>> Made solutions for math domain error. >>> >>> Prevented to take log of negative values, and division by zero. >>> >>> This though slows the implementation down. >>> >>> Systemtest Relax_disp.test_baldwin_synthetic_full went from 6.x seconds to >>> 8-9.x seconds. >>> >>> sr #3154: (https://gna.org/support/?3154) Implementation of Baldwin (2014) >>> B14 model - 2-site exact solution model for all time scales. >>> >>> This follows the tutorial for adding relaxation dispersion models at: >>> http://wiki.nmr-relax.com/Tutorial_for_adding_relaxation_dispersion_models_to_relax#Debugging >>> >>> Modified: >>> trunk/lib/dispersion/b14.py >>> >>> Modified: trunk/lib/dispersion/b14.py >>> URL: >>> http://svn.gna.org/viewcvs/relax/trunk/lib/dispersion/b14.py?rev=23091&r1=23090&r2=23091&view=diff >>> ============================================================================== >>> --- trunk/lib/dispersion/b14.py (original) >>> +++ trunk/lib/dispersion/b14.py Thu May 8 18:35:51 2014 >>> @@ -111,7 +111,7 @@ >>> >>> # Python module imports. >>> import numpy >>> -from numpy import arccosh, cos, cosh, log, sin, sinh, sqrt, power >>> +from numpy import arccosh, array, cos, cosh, in1d, log, nonzero, sin, >>> sinh, sqrt, power >>> >>> # Repetitive calculations (to speed up calculations). >>> g_fact = 1/sqrt(2) >>> @@ -201,43 +201,85 @@ >>> # Mixed term (complex) (E0 - iE2)/2. >>> E1 = (g3 - g4*1j) * tcp >>> >>> + # Complex. >>> + v1s = F0 * sinh(E0) - F2 * sin(E2)*1j >>> + >>> + # -2 * oG * t2. >>> + v4 = F1b * (-alpha_m - g3 ) + F1b * (dw - g4)*1j >>> + >>> + # Complex. >>> + ex1c = sinh(E1) >>> + >>> + # Off diagonal common factor. sinh fuctions. >>> + v5 = (-deltaR2 + kex + dw*1j) * v1s - 2. * (v4 + k_AB * F1a_plus_b) * >>> ex1c >>> + >>> # Real. The v_1c in paper. >>> v1c = F0 * cosh(E0) - F2 * cos(E2) >>> >>> - # Complex. >>> - v1s = F0 * sinh(E0) - F2 * sin(E2)*1j >>> - >>> - # Exact result for v2v3. >>> - v3 = sqrt(v1c**2 - 1.) >>> - >>> - # -2 * oG * t2. >>> - v4 = F1b * (-alpha_m - g3 ) + F1b * (dw - g4)*1j >>> - >>> - # Complex. >>> - ex1c = sinh(E1) >>> - >>> - # Off diagonal common factor. sinh fuctions. >>> - v5 = (-deltaR2 + kex + dw*1j) * v1s - 2. * (v4 + k_AB * F1a_plus_b) * >>> ex1c >>> - >>> - 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) ) >>> - >>> - # Loop over the time points, back calculating the R2eff values. >>> - for i in range(num_points): >>> - >>> - # The full formula. >>> - back_calc[i] = R2eff[i] >>> + # 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 >>> >>> >>> _______________________________________________ >>> relax (http://www.nmr-relax.com) >>> >>> This is the relax-commits mailing list >>> relax-comm...@gna.org >>> >>> 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@gna.org >> >> 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 relax-devel@gna.org 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
