Author: bugman
Date: Mon Jul 15 19:04:29 2013
New Revision: 20305
URL: http://svn.gna.org/viewcvs/relax?rev=20305&view=rev
Log:
Improvement of the error handling in the 'NS 2-site star' model.
The fA and pB parameters are no longer being checked. Instead a Mgx value of
0.0 is being checked
for. This catches additional problems. And now instead of the R2eff value
being set to zero, it is
set to 1e99. This is because log of zero is -inf, and then multiplied by a
negative constant gives
positive inf.
Modified:
branches/relax_disp/lib/dispersion/ns_2site_star.py
branches/relax_disp/target_functions/relax_disp.py
Modified: branches/relax_disp/lib/dispersion/ns_2site_star.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/ns_2site_star.py?rev=20305&r1=20304&r2=20305&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/ns_2site_star.py (original)
+++ branches/relax_disp/lib/dispersion/ns_2site_star.py Mon Jul 15 19:04:29
2013
@@ -39,7 +39,7 @@
from scipy.linalg import expm
-def r2eff_ns_2site_star(Rr=None, Rex=None, RCS=None, R=None, M0=None,
r20a=None, r20b=None, fA=None, pB=None, inv_tcpmg=None, tcp=None,
back_calc=None, num_points=None, power=None):
+def r2eff_ns_2site_star(Rr=None, Rex=None, RCS=None, R=None, M0=None,
r20a=None, r20b=None, fA=None, inv_tcpmg=None, tcp=None, back_calc=None,
num_points=None, power=None):
"""The 2-site numerical solution to the Bloch-McConnell equation using
complex conjugate matrices.
This function calculates and stores the R2eff values.
@@ -53,8 +53,6 @@
@type r20b: float
@keyword fA: The frequency of state A.
@type fA: float
- @keyword pB: The population of state B.
- @type pB: float
@keyword inv_tcpmg: The inverse of the total duration of the CPMG
element (in inverse seconds).
@type inv_tcpmg: float
@keyword tcp: The tau_CPMG times (1 / 4.nu1).
@@ -83,11 +81,6 @@
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
- # Catch zeros (to avoid matrix log failures).
- if fA == 0.0 and pB == 0.0:
- back_calc[i] = 0.0
- continue
-
# This matrix is a propagator that will evolve the magnetization
with the matrix R for a delay tcp.
expm_R_tcp = expm(R*tcp[i])
@@ -102,4 +95,7 @@
# The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
Mgx = Moft[0].real / M0[0]
- back_calc[i]= -inv_tcpmg * log(Mgx)
+ if Mgx == 0.0:
+ back_calc[i] = 1e99
+ else:
+ back_calc[i]= -inv_tcpmg * log(Mgx)
Modified: branches/relax_disp/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/target_functions/relax_disp.py?rev=20305&r1=20304&r2=20305&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Mon Jul 15 19:04:29
2013
@@ -554,15 +554,15 @@
dw_frq = dw[spin_index] * self.frqs[spin_index, frq_index]
# Back calculate the R2eff values.
- r2eff_ns_2site_star(Rr=self.Rr, Rex=self.Rex, RCS=self.RCS,
R=self.R, M0=self.M0, r20a=R20A[r20_index], r20b=R20B[r20_index], pB=pB,
fA=dw_frq, inv_tcpmg=self.inv_relax_time, tcp=self.tau_cpmg,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points, power=self.power)
-
- # For all missing data points, set the back-calculated value
to the measured values so that it has no effect on the chi-squared value.
- for point_index in range(self.num_disp_points):
- if self.missing[spin_index, frq_index, point_index]:
- self.back_calc[spin_index, frq_index, point_index] =
self.values[spin_index, frq_index, point_index]
-
- # Calculate and return the chi-squared value.
- chi2_sum += chi2(self.values[spin_index, frq_index],
self.back_calc[spin_index, frq_index], self.errors[spin_index, frq_index])
-
- # Return the total chi-squared value.
- return chi2_sum
+ r2eff_ns_2site_star(Rr=self.Rr, Rex=self.Rex, RCS=self.RCS,
R=self.R, M0=self.M0, r20a=R20A[r20_index], r20b=R20B[r20_index], fA=dw_frq,
inv_tcpmg=self.inv_relax_time, tcp=self.tau_cpmg,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points, power=self.power)
+
+ # For all missing data points, set the back-calculated value
to the measured values so that it has no effect on the chi-squared value.
+ for point_index in range(self.num_disp_points):
+ if self.missing[spin_index, frq_index, point_index]:
+ self.back_calc[spin_index, frq_index, point_index] =
self.values[spin_index, frq_index, point_index]
+
+ # Calculate and return the chi-squared value.
+ chi2_sum += chi2(self.values[spin_index, frq_index],
self.back_calc[spin_index, frq_index], self.errors[spin_index, frq_index])
+
+ # Return the total chi-squared value.
+ return chi2_sum
r20305 - in /branches/relax_disp: lib/dispersion/ns_2site_star.py target_functions/relax_disp.py