Author: bugman
Date: Fri Aug 30 11:45:03 2013
New Revision: 20725
URL: http://svn.gna.org/viewcvs/relax?rev=20725&view=rev
Log:
Created the 'NS R1rho 2-site' model target function.
This is the numerical model for the 2-site Bloch-McConnell equations for
R1rho data. The code
originates from the funNumrho.m file from the Skrynikov & Tollinger code (the
sim_all.tar file
https://gna.org/support/download.php?file_id=18404 attached to
https://gna.org/task/?7712#comment5).
This commit follows step 4 of the relaxation dispersion model addition
tutorial
(http://thread.gmane.org/gmane.science.nmr.relax.devel/3907).
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=20725&r1=20724&r2=20725&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Fri Aug 30 11:45:03
2013
@@ -37,12 +37,13 @@
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_r1rho_2site import ns_r1rho_2site
from lib.dispersion.ns_matrices import r180x_3d
from lib.dispersion.tp02 import r1rho_TP02
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 MODEL_CR72,
MODEL_CR72_FULL, MODEL_DPL94, MODEL_IT99, MODEL_LIST_FULL, MODEL_LM63,
MODEL_LM63_3SITE, MODEL_M61, MODEL_M61B, 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_R2EFF,
MODEL_TP02, MODEL_TSMFK01
+from specific_analyses.relax_disp.variables import MODEL_CR72,
MODEL_CR72_FULL, MODEL_DPL94, MODEL_IT99, MODEL_LIST_FULL, MODEL_LM63,
MODEL_LM63_3SITE, MODEL_M61, MODEL_M61B, 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:
@@ -189,6 +190,8 @@
if model in [MODEL_NS_CPMG_2SITE_3D, MODEL_NS_CPMG_2SITE_3D_FULL]:
self.M0 = zeros(7, float64)
self.M0[0] = 0.5
+ if model in [MODEL_NS_R1RHO_2SITE]:
+ self.M0 = zeros(6, float64)
# Some other data structures for the numerical solutions.
if model in [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]:
@@ -199,7 +202,8 @@
self.tau_cpmg[i] = 0.25 / self.cpmg_frqs[i]
self.power[i] = int(round(self.cpmg_frqs[i] *
self.relax_time))
- # The inverted relaxation delay.
+ # The inverted relaxation delay.
+ if model in [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
# Set up the model.
@@ -235,6 +239,8 @@
self.func = self.func_ns_cpmg_2site_star
if model == MODEL_TP02:
self.func = self.func_TP02
+ if model == MODEL_NS_R1RHO_2SITE:
+ self.func = self.func_ns_r1rho_2site
def calc_CR72_chi2(self, R20A=None, R20B=None, dw=None, pA=None,
kex=None):
@@ -929,6 +935,62 @@
return self.calc_ns_cpmg_2site_star_chi2(R20A=R20A, R20B=R20B,
dw=dw, pA=pA, kex=kex)
+ def func_ns_r1rho_2site(self, params):
+ """Target function for the reduced numerical solution for the 2-site
Bloch-McConnell equations for R1rho data.
+
+ @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.
+ r1rho_prime = params[:self.end_index[0]]
+ dw = params[self.end_index[0]:self.end_index[1]]
+ pA = params[self.end_index[1]]
+ kex = params[self.end_index[1]+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
+
+ # Chi-squared initialisation.
+ 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]
+
+ # Back calculate the R2eff values.
+ ns_r1rho_2site(M0=self.M0,
r1rho_prime=r1rho_prime[r20_index], omega=self.chemical_shifts[spin_index,
frq_index], r1=self.r1[spin_index, frq_index], pA=pA, pB=pB, dw=dw_frq,
k_AB=k_AB, k_BA=k_BA, spin_lock_fields=self.spin_lock_nu1,
relax_time=self.relax_time, inv_relax_time=self.inv_relax_time,
back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
+
+ # 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_TP02(self, params):
"""Target function for the Trott and Palmer (2002) R1rho
off-resonance 2-site model.
r20725 - /branches/relax_disp/target_functions/relax_disp.py