Author: bugman
Date: Fri Nov 22 15:58:32 2013
New Revision: 21601
URL: http://svn.gna.org/viewcvs/relax?rev=21601&view=rev
Log:
Big bug fixes for the 'MMQ 2-site' dispersion model.
The equations in Korzhnev et al., 2004 are incorrect! This is clear when
comparing the equations in
the paper to those implemented in the cpmg_fit source code. The equations in
relax have been
changed to match the correct one.
In addition, the matrix power factor must be found with the python
math.floor() function and not
int() as the later will sometimes round up.
Modified:
branches/relax_disp/lib/dispersion/mmq_2site.py
branches/relax_disp/target_functions/relax_disp.py
Modified: branches/relax_disp/lib/dispersion/mmq_2site.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/mmq_2site.py?rev=21601&r1=21600&r2=21601&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/mmq_2site.py (original)
+++ branches/relax_disp/lib/dispersion/mmq_2site.py Fri Nov 22 15:58:32 2013
@@ -30,6 +30,7 @@
"""
# Python module imports.
+from math import floor
from numpy import array, conj, dot, float64, log
# relax module imports.
@@ -62,7 +63,7 @@
matrix[1, 1] = -k_BA + 1.j*dw - R20B
-def r2eff_mmq_2site_mq(M0=None, F_vector=array([1, 0], float64), m1=None,
m2=None, R20A=None, R20B=None, pA=None, pB=None, dw=None, dwH=None,
k_AB=None, k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None,
num_points=None, power=None, n=None):
+def r2eff_mmq_2site_mq(M0=None, F_vector=array([1, 0], float64), m1=None,
m2=None, R20A=None, R20B=None, pA=None, pB=None, dw=None, dwH=None,
k_AB=None, k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None,
num_points=None, power=None):
"""The 2-site numerical solution to the Bloch-McConnell equation for MQ
data.
The notation used here comes from:
@@ -110,13 +111,11 @@
@type num_points: int
@keyword power: The matrix exponential power array.
@type power: numpy int16, rank-1 array
- @keyword n: The n value whereby one CPMG block is defined at
2n.
- @type n: numpy int16, rank-1 array
"""
# Populate the m1 and m2 matrices (only once per function call for
speed).
- populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, dw=dw+dwH, k_AB=k_AB,
k_BA=k_BA) # D+ matrix component.
- populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, dw=-dw+dwH, k_AB=k_AB,
k_BA=k_BA) # Z- matrix component.
+ populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, dw=-dw-dwH, k_AB=k_AB,
k_BA=k_BA) # D+ matrix component.
+ populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, dw=dw-dwH, k_AB=k_AB,
k_BA=k_BA) # Z- matrix component.
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
@@ -138,10 +137,24 @@
M1_M2_M2_M1_star = dot(M1_M2_star, M2_M1_star)
M2_M1_M1_M2_star = dot(M2_M1_star, M1_M2_star)
- # Matrices for even n.
- if n[i] % 2 == 0:
+ # Special case of 1 CPMG block - the power is zero.
+ if power[i] == 1:
+ # M2.M1
+ A = M2_M1
+
+ # M2*.M1*
+ B = M2_M1_star
+
+ # M1.M2
+ C = M1_M2
+
+ # M1*.M2*
+ D = M1_M2_star
+
+ # Matrices for even number of CPMG blocks.
+ elif power[i] % 2 == 0:
# The power factor (only calculate once).
- fact = int(n[i] / 2)
+ fact = int(floor(power[i] / 2))
# (M1.M2.M2.M1)^(n/2)
A = square_matrix_power(M1_M2_M2_M1, fact)
@@ -155,31 +168,31 @@
# (M1*.M2*.M2*.M1*)^(n/2)
D = square_matrix_power(M1_M2_M2_M1_star, fact)
- # Matrices for odd n.
+ # Matrices for odd number of CPMG blocks.
else:
# The power factor (only calculate once).
- fact = int((n[i] - 1) / 2)
-
- # (M1.M2.M2.M1)^((n-1)/2).M1.M2
+ fact = int(floor((power[i] - 1) / 2))
+
+ # M2.M1.(M1.M2.M2.M1)^((n-1)/2)
A = square_matrix_power(M1_M2_M2_M1, fact)
- A = dot(A, M1_M2)
-
- # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*
+ A = dot(M2_M1, A)
+
+ # M2*.M1*.(M1*.M2*.M2*.M1*)^((n-1)/2)
B = square_matrix_power(M1_M2_M2_M1_star, fact)
- B = dot(B, M1_M2_star)
-
- # (M2.M1.M1.M2)^((n-1)/2).M2.M1
+ B = dot(M2_M1_star, B)
+
+ # M1.M2.(M2.M1.M1.M2)^((n-1)/2)
C = square_matrix_power(M2_M1_M1_M2, fact)
- C = dot(C, M2_M1)
-
- # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*
+ C = dot(M1_M2, C)
+
+ # M1*.M2*.(M2*.M1*.M1*.M2*)^((n-1)/2)
D = square_matrix_power(M2_M1_M1_M2_star, fact)
- D = dot(D, M2_M1_star)
+ D = dot(M1_M2_star, D)
# The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
- A_B = dot(A, B)
- C_D = dot(C, D)
- Mx = dot(dot(F_vector, (A_B + C_D)), M0)
+ B_A = dot(B, A)
+ D_C = dot(D, C)
+ Mx = dot(dot(F_vector, (B_A + D_C)), M0)
Mx = Mx.real / 2.0
if Mx <= 0.0 or isNaN(Mx):
back_calc[i] = 1e99
@@ -187,7 +200,7 @@
back_calc[i]= -inv_tcpmg * log(Mx / pA)
-def r2eff_mmq_2site_sq_dq_zq(M0=None, F_vector=array([1, 0], float64),
m1=None, m2=None, R20A=None, R20B=None, pA=None, pB=None, dw=None, dwH=None,
k_AB=None, k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None,
num_points=None, power=None, n=None):
+def r2eff_mmq_2site_sq_dq_zq(M0=None, F_vector=array([1, 0], float64),
m1=None, m2=None, R20A=None, R20B=None, pA=None, pB=None, dw=None, dwH=None,
k_AB=None, k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None,
num_points=None, power=None):
"""The 2-site numerical solution to the Bloch-McConnell equation for SQ,
ZQ, and DQ data.
The notation used here comes from:
@@ -231,8 +244,6 @@
@type num_points: int
@keyword power: The matrix exponential power array.
@type power: numpy int16, rank-1 array
- @keyword n: The n value whereby one CPMG block is defined at
2n.
- @type n: numpy int16, rank-1 array
"""
# Populate the m1 and m2 matrices (only once per function call for
speed).
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=21601&r1=21600&r2=21601&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Fri Nov 22 15:58:32
2013
@@ -230,14 +230,6 @@
self.power[exp_type_index].append(zeros(self.num_disp_points[exp_type_index][frq_index],
int16))
for i in
range(self.num_disp_points[exp_type_index][frq_index]):
self.power[exp_type_index][frq_index][i] =
int(round(self.cpmg_frqs[exp_type_index][frq_index][i] *
self.relax_times[exp_type_index][frq_index]))
-
- # The n value.
- if model == MODEL_MMQ_2SITE:
- self.n = []
- for exp_type_index in range(len(values)):
- self.n.append([])
- for frq_index in range(len(values[exp_type_index][0])):
- self.n[exp_type_index].append(2 *
self.power[exp_type_index][frq_index])
# The tau_cpmg times - recalculated to avoid any user induced
truncation in the input files.
if model in [MODEL_MQ_CR72, MODEL_MMQ_2SITE, MODEL_NS_CPMG_2SITE_3D,
MODEL_NS_CPMG_2SITE_3D_FULL, MODEL_NS_CPMG_2SITE_EXPANDED,
MODEL_NS_CPMG_2SITE_STAR, MODEL_NS_CPMG_2SITE_STAR_FULL, MODEL_TSMFK01]:
@@ -968,7 +960,7 @@
aliased_dwH = dw_frq
# Back calculate the R2eff values for each experiment
type.
- self.r2eff_mmq[exp_index](M0=self.M0, m1=self.m1,
m2=self.m2, R20A=R20[r20_index], R20B=R20[r20_index], pA=pA, pB=pB,
dw=aliased_dw, dwH=aliased_dwH, k_AB=k_AB, k_BA=k_BA,
inv_tcpmg=self.inv_relax_times[exp_index][frq_index],
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index], n=self.n[exp_index][frq_index])
+ self.r2eff_mmq[exp_index](M0=self.M0, m1=self.m1,
m2=self.m2, R20A=R20[r20_index], R20B=R20[r20_index], pA=pA, pB=pB,
dw=aliased_dw, dwH=aliased_dwH, k_AB=k_AB, k_BA=k_BA,
inv_tcpmg=self.inv_relax_times[exp_index][frq_index],
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
# 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[exp_index][frq_index]):
r21601 - in /branches/relax_disp: lib/dispersion/mmq_2site.py target_functions/relax_disp.py