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 <[email protected]>: > 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, <[email protected]> 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 >> [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
