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
r23091 - /trunk/lib/dispersion/b14.py