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-commits@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