Author: bugman
Date: Tue Oct 15 18:54:44 2013
New Revision: 21126
URL: http://svn.gna.org/viewcvs/relax?rev=21126&view=rev
Log:
Created the 'MQ NS CPMG 2-site' model target function.
This is the Carver and Richards (1972) 2-site model expanded for MQ CPMG data
by Korzhnev et al.,
2004.
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:
branches/relax_disp/target_functions/relax_disp.py
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=21126&r1=21125&r2=21126&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Tue Oct 15 18:54:44
2013
@@ -35,6 +35,7 @@
from lib.dispersion.lm63_3site import r2eff_LM63_3site
from lib.dispersion.m61 import r1rho_M61
from lib.dispersion.m61b import r1rho_M61b
+from lib.dispersion.mq_cr72 import r2eff_mq_cr72
from lib.dispersion.mq_ns_cpmg_2site import r2eff_mq_ns_cpmg_2site
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
@@ -45,7 +46,7 @@
from lib.dispersion.tsmfk01 import r2eff_TSMFK01
from lib.errors import RelaxError
from target_functions.chi2 import chi2
-from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG,
EXP_TYPE_MQ_CPMG, EXP_TYPE_MQ_R1RHO, EXP_TYPE_R1RHO, MODEL_CR72,
MODEL_CR72_FULL, MODEL_DPL94, MODEL_IT99, MODEL_LIST_CPMG, MODEL_LIST_FULL,
MODEL_LIST_MQ_CPMG, MODEL_LIST_MQ_R1RHO, MODEL_LIST_R1RHO, MODEL_LM63,
MODEL_LM63_3SITE, MODEL_M61, MODEL_M61B, MODEL_MQ_NS_CPMG_2SITE, 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_TP02,
MODEL_TSMFK01
+from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG,
EXP_TYPE_MQ_CPMG, EXP_TYPE_MQ_R1RHO, EXP_TYPE_R1RHO, MODEL_CR72,
MODEL_CR72_FULL, MODEL_DPL94, MODEL_IT99, MODEL_LIST_CPMG, MODEL_LIST_FULL,
MODEL_LIST_MQ_CPMG, MODEL_LIST_MQ_R1RHO, MODEL_LIST_R1RHO, MODEL_LM63,
MODEL_LM63_3SITE, MODEL_M61, MODEL_M61B, MODEL_MQ_CR72,
MODEL_MQ_NS_CPMG_2SITE, 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_TP02, MODEL_TSMFK01
class Dispersion:
@@ -173,7 +174,7 @@
# The spin and dependent parameters (phi_ex, dw, padw2).
self.end_index.append(self.end_index[-1] + self.num_spins)
- if model in [MODEL_IT99, MODEL_LM63_3SITE, MODEL_MQ_NS_CPMG_2SITE]:
+ if model in [MODEL_IT99, MODEL_LM63_3SITE, MODEL_MQ_CR72,
MODEL_MQ_NS_CPMG_2SITE]:
self.end_index.append(self.end_index[-1] + self.num_spins)
# Set up the matrices for the numerical solutions.
@@ -204,14 +205,14 @@
self.M0 = zeros(6, float64)
# Some other data structures for the analytical and numerical
solutions.
- if model in [MODEL_MQ_NS_CPMG_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_MQ_NS_CPMG_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]:
# The tau_cpmg times.
self.tau_cpmg = zeros(self.num_disp_points, float64)
for i in range(self.num_disp_points):
self.tau_cpmg[i] = 0.25 / self.cpmg_frqs[i]
# Some other data structures for the numerical solutions.
- if model in [MODEL_MQ_NS_CPMG_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]:
+ if model in [MODEL_MQ_CR72, MODEL_MQ_NS_CPMG_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]:
# The matrix exponential power array.
self.power = zeros(self.num_disp_points, int16)
for i in range(self.num_disp_points):
@@ -223,7 +224,7 @@
self.spin_lock_omega1_squared = self.spin_lock_omega1 ** 2
# The inverted relaxation delay.
- if model in [MODEL_MQ_NS_CPMG_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_MQ_NS_CPMG_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]:
self.inv_relax_time = 1.0 / relax_time
# Special matrices for the multi-quantum CPMG 2-site numerical model.
@@ -266,6 +267,8 @@
self.func = self.func_TP02
if model == MODEL_NS_R1RHO_2SITE:
self.func = self.func_ns_r1rho_2site
+ if model == MODEL_MQ_CR72:
+ self.func = self.func_mq_CR72
if model == MODEL_MQ_NS_CPMG_2SITE:
self.func = self.func_mq_ns_cpmg_2site
@@ -795,8 +798,8 @@
return chi2_sum
- def func_mq_ns_cpmg_2site(self, params):
- """Target function for the Ishima and Torchia (1999) 2-site model
for all timescales with pA >> pB.
+ def func_mq_CR72(self, params):
+ """Target function for the CR72 model extended for MQ CPMG data.
@param params: The vector of parameter values.
@type params: numpy rank-1 float array
@@ -839,6 +842,64 @@
dwH_frq = dwH[spin_index] * self.frqs[spin_index, frq_index]
# Back calculate the R2eff values.
+ r2eff_mq_cr72(r20=R20[r20_index], pA=pA, pB=pB, dw=dw_frq,
dwH=dwH_frq, kex=kex, k_AB=k_AB, k_BA=k_BA, 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
+
+
+ def func_mq_ns_cpmg_2site(self, params):
+ """Target function for the Ishima and Torchia (1999) 2-site model
for all timescales with pA >> pB.
+
+ @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 = params[self.end_index[0]:self.end_index[1]]
+ dwH = params[self.end_index[1]:self.end_index[2]]
+ pA = params[self.end_index[2]]
+ kex = params[self.end_index[2]+1]
+
+ # Once off parameter conversions.
+ pB = 1.0 - pA
+ k_BA = pA * kex
+ k_AB = pB * kex
+
+ # 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
+
+ # Initialise.
+ chi2_sum = 0.0
+
+ # 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 + spin_index*self.num_frq
+
+ # Convert dw from ppm to rad/s.
+ dw_frq = dw[spin_index] * self.frqs[spin_index, frq_index]
+ dwH_frq = dwH[spin_index] * self.frqs[spin_index, frq_index]
+
+ # Back calculate the R2eff values.
r2eff_mq_ns_cpmg_2site(M0=self.M0, m1=self.m1, m2=self.m2,
r20=R20[r20_index], pA=pA, pB=pB, dw=dw_frq, dwH=dwH_frq, k_AB=k_AB,
k_BA=k_BA, 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.
r21126 - /branches/relax_disp/target_functions/relax_disp.py