Author: bugman
Date: Tue Dec 3 17:23:43 2013
New Revision: 21753
URL: http://svn.gna.org/viewcvs/relax?rev=21753&view=rev
Log:
Created the target functions for the 'NS MMQ 3-site' models.
This is for the 'NS MMQ 3-site' and 'NS MMQ 3-site (linear)' CPMG dispersion
models.
This follows the tutorial for adding relaxation dispersion models at:
http://wiki.nmr-relax.com/Tutorial_for_adding_relaxation_dispersion_models_to_relax#The_target_function
Modified:
trunk/target_functions/relax_disp.py
Modified: trunk/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/relax_disp.py?rev=21753&r1=21752&r2=21753&view=diff
==============================================================================
--- trunk/target_functions/relax_disp.py (original)
+++ trunk/target_functions/relax_disp.py Tue Dec 3 17:23:43 2013
@@ -42,6 +42,7 @@
from lib.dispersion.ns_cpmg_2site_3d import r2eff_ns_cpmg_2site_3D
from lib.dispersion.ns_cpmg_2site_expanded import
r2eff_ns_cpmg_2site_expanded
from lib.dispersion.ns_cpmg_2site_star import r2eff_ns_cpmg_2site_star
+from lib.dispersion.ns_mmq_3site import r2eff_ns_mmq_3site_mq,
r2eff_ns_mmq_3site_sq_dq_zq
from lib.dispersion.ns_r1rho_2site import ns_r1rho_2site
from lib.dispersion.ns_matrices import r180x_3d
from lib.dispersion.tp02 import r1rho_TP02
@@ -50,7 +51,7 @@
from lib.errors import RelaxError
from lib.float import isNaN
from target_functions.chi2 import chi2
-from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG_DQ,
EXP_TYPE_CPMG_MQ, EXP_TYPE_CPMG_PROTON_MQ, EXP_TYPE_CPMG_PROTON_SQ,
EXP_TYPE_CPMG_SQ, EXP_TYPE_CPMG_ZQ, EXP_TYPE_R1RHO, MODEL_CR72,
MODEL_CR72_FULL, MODEL_DPL94, MODEL_IT99, MODEL_LIST_CPMG, MODEL_LIST_FULL,
MODEL_LIST_MMQ, MODEL_LIST_MQ_CPMG, MODEL_LIST_R1RHO, MODEL_LM63,
MODEL_LM63_3SITE, MODEL_M61, MODEL_M61B, MODEL_MMQ_2SITE, MODEL_MP05,
MODEL_MQ_CR72, MODEL_NOREX, 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_NS_R1RHO_2SITE, MODEL_R2EFF, MODEL_TAP03, MODEL_TP02, MODEL_TSMFK01
+from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG_DQ,
EXP_TYPE_CPMG_MQ, EXP_TYPE_CPMG_PROTON_MQ, EXP_TYPE_CPMG_PROTON_SQ,
EXP_TYPE_CPMG_SQ, EXP_TYPE_CPMG_ZQ, EXP_TYPE_R1RHO, MODEL_CR72,
MODEL_CR72_FULL, MODEL_DPL94, MODEL_IT99, MODEL_LIST_CPMG, MODEL_LIST_FULL,
MODEL_LIST_MMQ, MODEL_LIST_MQ_CPMG, MODEL_LIST_R1RHO, MODEL_LM63,
MODEL_LM63_3SITE, MODEL_M61, MODEL_M61B, MODEL_MMQ_2SITE, MODEL_MP05,
MODEL_MQ_CR72, MODEL_NOREX, 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_NS_MMQ_3SITE,
MODEL_NS_MMQ_3SITE_LINEAR, MODEL_NS_R1RHO_2SITE, MODEL_R2EFF, MODEL_TAP03,
MODEL_TP02, MODEL_TSMFK01
class Dispersion:
@@ -84,6 +85,8 @@
- 'NS CPMG 2-site star full': The full numerical solution for
the 2-site Bloch-McConnell equations for CPMG data using complex conjugate
matrices.
- 'NS CPMG 2-site expanded': The numerical solution for the
2-site Bloch-McConnell equations for CPMG data expanded using Maple by
Nikolai Skrynnikov.
- 'MMQ 2-site': The numerical solution for the 2-site
Bloch-McConnell equations for combined proton-heteronuclear SQ, ZD, DQ, and
MQ CPMG data with R20A = R20B.
+ - 'NS MMQ 3-site linear': The numerical solution for the 3-site
Bloch-McConnell equations linearised with kAC = kCA = 0 for combined
proton-heteronuclear SQ, ZD, DQ, and MQ CPMG data with R20A = R20B = R20C.
+ - 'NS MMQ 3-site': The numerical solution for the 3-site
Bloch-McConnell equations for combined proton-heteronuclear SQ, ZD, DQ, and
MQ CPMG data with R20A = R20B = R20C.
@keyword model: The relaxation dispersion model to fit.
@@ -195,6 +198,10 @@
self.end_index.append(self.end_index[-1] + self.num_spins)
if model in [MODEL_IT99, MODEL_LM63_3SITE, MODEL_MQ_CR72,
MODEL_MMQ_2SITE]:
self.end_index.append(self.end_index[-1] + self.num_spins)
+ elif model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR]:
+ self.end_index.append(self.end_index[-1] + self.num_spins)
+ self.end_index.append(self.end_index[-1] + self.num_spins)
+ self.end_index.append(self.end_index[-1] + self.num_spins)
# Set up the matrices for the numerical solutions.
if model in [MODEL_NS_CPMG_2SITE_STAR,
MODEL_NS_CPMG_2SITE_STAR_FULL]:
@@ -217,6 +224,8 @@
# This is a vector that contains the initial magnetizations
corresponding to the A and B state transverse magnetizations.
if model in [MODEL_MMQ_2SITE, MODEL_NS_CPMG_2SITE_STAR,
MODEL_NS_CPMG_2SITE_STAR_FULL]:
self.M0 = zeros(2, float64)
+ if model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR]:
+ self.M0 = zeros(3, float64)
if model in [MODEL_NS_CPMG_2SITE_3D, MODEL_NS_CPMG_2SITE_3D_FULL]:
self.M0 = zeros(7, float64)
self.M0[0] = 0.5
@@ -224,7 +233,7 @@
self.M0 = zeros(6, float64)
# Special CPMG-type data structures.
- 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]:
+ 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_NS_MMQ_3SITE,
MODEL_NS_MMQ_3SITE_LINEAR, MODEL_TSMFK01]:
# The number of CPMG blocks.
self.power = []
for exp_type_index in range(self.num_exp):
@@ -267,13 +276,16 @@
self.spin_lock_omega1_squared[exp_type_index][frq_index]
= self.spin_lock_omega1[exp_type_index][frq_index] ** 2
# The inverted relaxation delay.
- 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_NS_R1RHO_2SITE]:
+ 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_NS_MMQ_3SITE,
MODEL_NS_MMQ_3SITE_LINEAR, MODEL_NS_R1RHO_2SITE]:
self.inv_relax_times = 1.0 / relax_times
- # Special storage matrices for the multi-quantum CPMG 2-site
numerical model.
+ # Special storage matrices for the multi-quantum CPMG N-site
numerical models.
if model == MODEL_MMQ_2SITE:
self.m1 = zeros((2, 2), complex64)
self.m2 = zeros((2, 2), complex64)
+ elif model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR]:
+ self.m1 = zeros((3, 3), complex64)
+ self.m2 = zeros((3, 3), complex64)
# Set up the model.
if model == MODEL_NOREX:
@@ -318,6 +330,10 @@
self.func = self.func_mq_CR72
if model == MODEL_MMQ_2SITE:
self.func = self.func_mmq_2site
+ if model == MODEL_NS_MMQ_3SITE:
+ self.func = self.func_ns_mmq_3site
+ if model == MODEL_NS_MMQ_3SITE_LINEAR:
+ self.func = self.func_ns_mmq_3site_linear
def calc_CR72_chi2(self, R20A=None, R20B=None, dw=None, pA=None,
kex=None):
@@ -479,6 +495,121 @@
return chi2_sum
+ def calc_ns_mmq_3site_chi2(self, R20A=None, R20B=None, R20C=None,
dw_AB=None, dw_BC=None, dw_AC=None, dwH_AB=None, dwH_BC=None, dwH_AC=None,
pA=None, pB=None, kex_AB=None, kex_BC=None, kex_AC=None):
+ """Calculate the chi-squared value for the 'NS MMQ 3-site' models.
+
+ @keyword R20A: The R2 value for state A in the absence of
exchange.
+ @type R20A: list of float
+ @keyword R20B: The R2 value for state B in the absence of
exchange.
+ @type R20B: list of float
+ @keyword R20C: The R2 value for state C in the absence of
exchange.
+ @type R20C: list of float
+ @keyword dw_AB: The chemical exchange difference between states
A and B in rad/s.
+ @type dw_AB: float
+ @keyword dw_BC: The chemical exchange difference between states
B and C in rad/s.
+ @type dw_BC: float
+ @keyword dwH_AB: The proton chemical exchange difference between
states A and B in rad/s.
+ @type dwH_AB: float
+ @keyword dwH_BC: The proton chemical exchange difference between
states B and C in rad/s.
+ @type dwH_BC: float
+ @keyword pA: The population of state A.
+ @type pA: float
+ @keyword kex_AB: The rate of exchange between states A and B.
+ @type kex_AB: float
+ @keyword kex_BC: The rate of exchange between states B and C.
+ @type kex_BC: float
+ @keyword kex_AC: The rate of exchange between states A and C.
+ @type kex_AC: float
+ @return: The chi-squared value.
+ @rtype: float
+ """
+
+ # Once off parameter conversions.
+ pC = 1.0 - pA - pB
+ k_BA = pA * kex_AB
+ k_AB = pB * kex_AB
+ k_CA = pA * kex_AC
+ k_AC = pC * kex_AC
+ k_CB = pB * kex_BC
+ k_BC = pC * kex_BC
+ dw_AC = dw_AB + dw_BC
+ dwH_AC = dwH_AB + dwH_BC
+
+ # This is a vector that contains the initial magnetizations
corresponding to the A and B state transverse magnetizations.
+ self.M0[0] = pA
+ self.M0[1] = pB
+ self.M0[2] = pC
+
+ # Initialise.
+ chi2_sum = 0.0
+
+ # Loop over the experiment types.
+ for exp_index in range(self.num_exp):
+ # Loop over the spins.
+ for spin_index in range(self.num_spins):
+ # Loop over the spectrometer frequencies.
+ for frq_index in range(self.num_frq):
+ # The R20 index.
+ r20_index = frq_index + exp_index*self.num_frq +
spin_index*self.num_frq*self.num_exp
+
+ # Convert dw from ppm to rad/s.
+ dw_AB_frq = dw_AB[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
+ dw_BC_frq = dw_BC[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
+ dw_AC_frq = dw_AC[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
+ dwH_AB_frq = dwH_AB[spin_index] *
self.frqs_H[exp_index][spin_index][frq_index]
+ dwH_BC_frq = dwH_BC[spin_index] *
self.frqs_H[exp_index][spin_index][frq_index]
+ dwH_AC_frq = dwH_AC[spin_index] *
self.frqs_H[exp_index][spin_index][frq_index]
+
+ # Alias the dw frequency combinations.
+ aliased_dwH_AB = 0.0
+ aliased_dwH_BC = 0.0
+ aliased_dwH_AC = 0.0
+ if self.exp_types[exp_index] == EXP_TYPE_CPMG_SQ:
+ aliased_dw_AB = dw_AB_frq
+ aliased_dw_BC = dw_BC_frq
+ aliased_dw_AC = dw_AC_frq
+ elif self.exp_types[exp_index] ==
EXP_TYPE_CPMG_PROTON_SQ:
+ aliased_dw_AB = dwH_AB_frq
+ aliased_dw_BC = dwH_BC_frq
+ aliased_dw_AC = dwH_AC_frq
+ elif self.exp_types[exp_index] == EXP_TYPE_CPMG_DQ:
+ aliased_dw_AB = dw_AB_frq + dwH_AB_frq
+ aliased_dw_BC = dw_BC_frq + dwH_BC_frq
+ aliased_dw_AC = dw_AC_frq + dwH_AC_frq
+ elif self.exp_types[exp_index] == EXP_TYPE_CPMG_ZQ:
+ aliased_dw_AB = dw_AB_frq - dwH_AB_frq
+ aliased_dw_BC = dw_BC_frq - dwH_BC_frq
+ aliased_dw_AC = dw_AC_frq - dwH_AC_frq
+ elif self.exp_types[exp_index] == EXP_TYPE_CPMG_MQ:
+ aliased_dw_AB = dw_AB_frq
+ aliased_dw_BC = dw_BC_frq
+ aliased_dw_AC = dw_AC_frq
+ aliased_dwH_AB = dwH_AB_frq
+ aliased_dwH_BC = dwH_BC_frq
+ aliased_dwH_AC = dwH_AC_frq
+ elif self.exp_types[exp_index] ==
EXP_TYPE_CPMG_PROTON_MQ:
+ aliased_dw_AB = dwH_AB_frq
+ aliased_dw_BC = dwH_BC_frq
+ aliased_dw_AC = dwH_AC_frq
+ aliased_dwH_AB = dw_AB_frq
+ aliased_dwH_BC = dw_BC_frq
+ aliased_dwH_AC = dw_AC_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=R20A[r20_index], R20B=R20B[r20_index], R20C=R20C[r20_index],
pA=pA, pB=pB, pC=pC, dw_AB=aliased_dw_AB, dw_BC=aliased_dw_BC,
dw_AC=aliased_dw_AC, dwH_AB=aliased_dwH_AB, dwH_BC=aliased_dwH_BC,
dwH_AC=aliased_dwH_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC,
k_CA=k_CA, 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]):
+ if
self.missing[exp_index][spin_index][frq_index][point_index]:
+
self.back_calc[exp_index][spin_index][frq_index][point_index] =
self.values[exp_index][spin_index][frq_index][point_index]
+
+ # Calculate and return the chi-squared value.
+ chi2_sum +=
chi2(self.values[exp_index][spin_index][frq_index],
self.back_calc[exp_index][spin_index][frq_index],
self.errors[exp_index][spin_index][frq_index])
+
+ # Return the total chi-squared value.
+ return chi2_sum
+
+
def experiment_type_setup(self):
"""Check the experiment types and simplify data structures.
@@ -492,11 +623,22 @@
if self.model in MODEL_LIST_MMQ:
# Alias the r2eff functions.
self.r2eff_mmq = []
+
+ # Loop over the experiment types.
for exp_index in range(self.num_exp):
+ # SQ, DQ and ZQ data types.
if self.exp_types[exp_index] in [EXP_TYPE_CPMG_SQ,
EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_DQ, EXP_TYPE_CPMG_ZQ]:
- self.r2eff_mmq.append(r2eff_mmq_2site_sq_dq_zq)
+ if self.model == MODEL_MMQ_2SITE:
+ self.r2eff_mmq.append(r2eff_mmq_2site_sq_dq_zq)
+ else:
+ self.r2eff_mmq.append(r2eff_ns_mmq_3site_sq_dq_zq)
+
+ # MQ data types.
elif self.exp_types[exp_index] in [EXP_TYPE_CPMG_MQ,
EXP_TYPE_CPMG_PROTON_MQ]:
- self.r2eff_mmq.append(r2eff_mmq_2site_mq)
+ if self.model == MODEL_MMQ_2SITE:
+ self.r2eff_mmq.append(r2eff_mmq_2site_mq)
+ else:
+ self.r2eff_mmq.append(r2eff_ns_mmq_3site_mq)
# The single data type models.
else:
@@ -1250,6 +1392,63 @@
return self.calc_ns_cpmg_2site_star_chi2(R20A=R20A, R20B=R20B,
dw=dw, pA=pA, kex=kex)
+ def func_ns_mmq_3site(self, params):
+ """Target function for the combined SQ, ZQ, DQ and MQ 3-site MMQ
CPMG numeric solution.
+
+ @param params: The vector of parameter values.
+ @type params: numpy rank-1 float array
+ @return: The chi-squared value.
+ @rtype: float
+ """
+
+ # Scaling.
+ if self.scaling_flag:
+ params = dot(params, self.scaling_matrix)
+
+ # Unpack the parameter values.
+ R20 = params[:self.end_index[0]]
+ dw_AB = params[self.end_index[0]:self.end_index[1]]
+ dw_BC = params[self.end_index[1]:self.end_index[2]]
+ dwH_AB = params[self.end_index[2]:self.end_index[3]]
+ dwH_BC = params[self.end_index[3]:self.end_index[4]]
+ pA = params[self.end_index[4]]
+ kex_AB = params[self.end_index[4]+1]
+ pB = params[self.end_index[4]+2]
+ kex_BC = params[self.end_index[4]+3]
+ kex_AC = params[self.end_index[4]+4]
+
+ # Calculate and return the chi-squared value.
+ return self.calc_ns_mmq_3site_chi2(R20A=R20, R20B=R20, R20C=R20,
dw_AB=dw_AB, dw_BC=dw_BC, dwH_AB=dwH_AB, dwH_BC=dwH_BC, pA=pA, pB=pB,
kex_AB=kex_AB, kex_BC=kex_BC, kex_AC=kex_AC)
+
+
+ def func_ns_mmq_3site_linear(self, params):
+ """Target function for the combined SQ, ZQ, DQ and MQ 3-site
linearised MMQ CPMG numeric solution.
+
+ @param params: The vector of parameter values.
+ @type params: numpy rank-1 float array
+ @return: The chi-squared value.
+ @rtype: float
+ """
+
+ # Scaling.
+ if self.scaling_flag:
+ params = dot(params, self.scaling_matrix)
+
+ # Unpack the parameter values.
+ R20 = params[:self.end_index[0]]
+ dw_AB = params[self.end_index[0]:self.end_index[1]]
+ dw_BC = params[self.end_index[1]:self.end_index[2]]
+ dwH_AB = params[self.end_index[2]:self.end_index[3]]
+ dwH_BC = params[self.end_index[3]:self.end_index[4]]
+ pA = params[self.end_index[4]]
+ kex_AB = params[self.end_index[4]+1]
+ pB = params[self.end_index[4]+2]
+ kex_BC = params[self.end_index[4]+3]
+
+ # Calculate and return the chi-squared value.
+ return self.calc_ns_mmq_3site_chi2(R20A=R20, R20B=R20, R20C=R20,
dw_AB=dw_AB, dw_BC=dw_BC, dwH_AB=dwH_AB, dwH_BC=dwH_BC, pA=pA, pB=pB,
kex_AB=kex_AB, kex_BC=kex_BC, kex_AC=0.0)
+
+
def func_ns_r1rho_2site(self, params):
"""Target function for the reduced numerical solution for the 2-site
Bloch-McConnell equations for R1rho data.
r21753 - /trunk/target_functions/relax_disp.py