r20492 - /branches/relax_disp/target_functions/relax_disp.py -- July 31, 2013 - 16:56

Author: bugman
Date: Wed Jul 31 16:56:39 2013
New Revision: 20492

URL: http://svn.gna.org/viewcvs/relax?rev=20492&view=rev
Log:
Created the 'TP02' model target function.

This is the Trott and Palmer 2002 R1rho analytic model for 2-site exchange.

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=20492&r1=20491&r2=20492&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Wed Jul 31 16:56:39 
2013
@@ -38,9 +38,10 @@
 from lib.dispersion.ns_2site_expanded import r2eff_ns_2site_expanded
 from lib.dispersion.ns_2site_star import r2eff_ns_2site_star
 from lib.dispersion.ns_matrices import r180x_3d
+from lib.dispersion.tp02 import r1rho_TP02
 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_2SITE_3D, 
MODEL_NS_2SITE_3D_FULL, MODEL_NS_2SITE_EXPANDED, MODEL_NS_2SITE_STAR, 
MODEL_NS_2SITE_STAR_FULL, MODEL_R2EFF
+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_2SITE_3D, 
MODEL_NS_2SITE_3D_FULL, MODEL_NS_2SITE_EXPANDED, MODEL_NS_2SITE_STAR, 
MODEL_NS_2SITE_STAR_FULL, MODEL_R2EFF, MODEL_TP02
 
 
 class Dispersion:
@@ -61,6 +62,7 @@
             - 'M61':  The Meiboom (1961) 2-site fast exchange model for 
R1rho-type experiments.
             - 'DPL94':  The Davis, Perlman and London (1994) 2-site fast 
exchange model for R1rho-type experiments.
             - 'M61 skew':  The Meiboom (1961) on-resonance 2-site model with 
skewed populations (pA >> pB) for R1rho-type experiments.
+            - 'TP02':  The Trott and Palmer (2002) 2-site exchange model for 
R1rho-type experiments.
 
         The following numerical models are currently supported:
 
@@ -210,6 +212,8 @@
             self.func = self.func_ns_2site_star_full
         if model == MODEL_NS_2SITE_STAR:
             self.func = self.func_ns_2site_star
+        if model == MODEL_TP02:
+            self.func = self.func_TP02
 
 
     def calc_CR72_chi2(self, R20A=None, R20B=None, dw=None, pA=None, 
kex=None):
@@ -902,3 +906,53 @@
 
         # Calculate and return the chi-squared value.
         return self.calc_ns_2site_star_chi2(R20A=R20A, R20B=R20B, dw=dw, 
pA=pA, kex=kex)
+
+
+    def func_TP02(self, params):
+        """Target function for the Trott and Palmer (2002) R1rho 
off-resonance 2-site model.
+
+        @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]]
+        pA = params[self.end_index[1]]
+        kex = params[self.end_index[1]+1]
+
+        # Once off parameter conversions.
+        pB = 1.0 - pA
+
+        # 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]
+
+                # Back calculate the R1rho values.
+                r1rho_TP02(r1rho_prime=R20[r20_index], pA=pA, pB=pB, 
dw=dw_frq, kex=kex, spin_lock_fields=self.spin_lock_nu1, 
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